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Abstract 



lis paper provides, from one side, a review of the theory of the cosmological mass function from a theoretical 
po^^ of view, starting from the seminal paper of Press & Shechter (1974) to the last developments (Del Popolo 
& g^mbera (1998, 1999), Sheth k Tormen 1999 (ST), Sheth, Mo k Tormen 2001 (STl), Jenkins et al. 2001 
(JOO), Shet k Tormen 2002 (ST2), Del Popolo 2002a, Yagi et al. 2004 (YNY)), and from another side some 
improvements on the multiplicity function models in literature. For what concerns this second aspect, I compare 
the^umerical multiplicity function given in Yahagi, Nagashima k Yoshii (2004), with the theoretical multiplicity 
function obtained in the present paper by means of the excursion set model and an improved version of the barrier 
sh^^ obtained in Del Popolo k Gambera (1998), which implicitly takes account of total angular momentum 
ac^^Qred by the proto-structure during evolution and of a non-zero cosmological constant. I show that the 
mijjbiplicity function obtained in the present paper, is in better agreement with Yahagi, Nagashima k Yoshii 
(2^^) simulations than other previous models (Sheth k Tormen 1999; Sheth, Mo k Tormen 2001; Sheth k 
To®ien 2002; Jenkins et al. 2001). The multiplicity function of the present paper gives a good fit to simulations 
rp stTi|-. s as the fit function proposed by Yahagi, Nagashima k Yoshii (2004), but differently from that it was 
obt^ifted from a sound theoretical background. Then, I calculate the mass function evolution in a ACDM model 
by ^eans of the previous model. I compare the result with Reed et al. (2003) (R03), who used a high resolution 
AC^M numerical simulation to calculate the mass function of dark matter haloes down to the scale of dwarf 
galaxies, back to a redshift of fifteen. I show that the mass function obtained in the present paper, gives similar 
pre^ctions to the Sheth k Tormen mass function but it does not show the overprediction of extremely rare 
bb^^ts shown by the Sheth and Tormen mass function. The results confirm previous findings that the simulated 
halSSmass function can be described solely by the variance of the mass distribution, and thus has no explicit 
redshift dependence. I moreover show that the PS-like approach together with the ellipsoidal model introduced 
in Del Popolo (2002b) gives a better description of the theoretical mass function. 

Subject headings: clusters of galaxies — cosmology: theory — dark matter — galaxies: clustering — large- 
scale structure of the universe 



1. Introduction 

The universe tliat we observe appears quite clumpy and inhomogeneous on a spatial scale of ~ 200h~^ Mpc. Beyond 
that scale mass clumps appear to be homogeneously distributed. We observe massive clumps such as galaxies, groups of 
galaxies, clusters of galaxies, super-clusters (the most massive one among the hierarchy of structures) fill the space spanning 
a wide range of mass scales. The respective mass ranges characterizing these systems, approximately, are ~ 10^ — lO^^M©, 
~ 1O"M0, ~ lO^^-lO^^'A/o, and ~ IO^^Mq where Mq is the solar mass. All these objects, combinedly, form the structure 
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what is known as the Large Scale Structure (LSS) in the universe. One of the most fundamental challenges in present 
universe is to understand the formation and evolution of the LSS. In order to understand the LSS, we need to have a 
theoretical framework within which predictions for structure formation can be made. 

The leading idea of all structure formation theories is that structures was born from small perturbations in the otherwise 
imiform distribution of matter in the early Universe, which is supposed to be, in great part, dark (matter not detectable 

through light emission). 

With the term Dark Matter cosmologists indicate an hypothetic material component of the universe which does not 
emit directly electromagnetic radiation (unless it decays in particles having this property (Sciama 1990, but also see 
Bowyer et al. 1999)). 

Dark matter, cannot be revealed directly, but nevertheless it is necessary to postulate its existence in order to explain 
the discrepancies between the observed dynamical properties of galaxies and clusters of galaxies and the theoretical 
predictions based upon models of these objects assuming that the only matter present is the visible one. The original 
hypotheses on Dark Matter go back to measures performed by Oort (1932) of the surface density of matter in the galactic 
disk, which was obtained through the study of the stars motion in direction orthogonal to the galactic plane. The result 
obtained by Oort, which was after him named "Oort Limit", gave a value of p = 0.15Mopc~^ for the mass density, and 
a mass, in the region studied, superior to that present in stars. Nowadays, we know that the quoted discrepancy is due 
to the presence of HI in the solar neighborhood. Other studies (Zwicky 1933; Smith 1936) showed the existence of a 
noteworthy discrepancy between the virial mass of clusters (e.g. Coma Cluster) and the total mass contained in galax;ies 
of the same clusters. These and other researches from the thirties to now, have confirmed that a great part of the mass 
in the universe does not emit radiation that can be directly observed. 

The study of Dark Matter has as its finality the explanation of formation of galaxies and in general of cosmic 
structures. For this reason, in the last decades, the origin of cosmic structures has been "framed" in models in which 
Dark Matter constitutes the skeleton of cosmic structures and supply the most part of the mass of which the same is made. 

The mass distribution of cosmic structures, i.e. the number of objects per unit volume and unit mass interval, is 
commonly called mass function or multiplicity function; it will be referred to as MF throughout the text. The determination 

of the cosmological MF is a difficult, not fully solved problem, both from the theoretical and the observational point of view. 
An analytical exact prediction of the MF is hampered, even in the simplest cosmological models, by the fact that highly 
non-linear gravitational dynamics is involved in the formation of high density objects; it is well known that the gravitational 
collapse problem has never been exactly solved, except in the case of simple symmetries (spherical planar). Large N-body 
simulations can be used to determine the MF of simulated halos; however such simulations, besides being time-expensive 
and limited in resolution, provide a numerical estimate of the final solution of the problem without directly shedding light 
on the difficult problem of gravitational collapse. Approximate analytical arguments, while being of limited validity, can 
provide useful and fully-understood solutions, which can then be compared to the results of N-body simulations. 

There is general consensus in setting the birth date of the MF theory in 1974, when the seminal paper of Press & 
Schechter (hereafter PS) was published (the same PS formula can be found in Doroshkevich 1967). That paper proposed a 
heuristic procedure, based on linear theory, to obtain the distribution of the masses of collapsed clumps. That work inspired 
the fit for the galaxy luminosity function proposed by Schechter (1976), but received a limited attention for more than a 
decade. A real explosion of attention to the MF theory started in 1988, when the first large N-body simulations started 
to reveal a surprising adherence of their results with the PS formula. Many authors tried to extend the PS procedure in 
many directions, or proposed different, alternative or complementary procedures (Bond, Cole, Efstathiou & Kaiser (1991) 
(BCEK)). These last years are witnessing a new wave of interest, which has not yet been exhausted (see for example ST, 
Jenkins et al. 2001 (JOl), STl, ST2, Del Popolo 2002a, Yahagi et al. 2004 (YNY), Del Popolo 2005). 

As mentioned before, the original PS work was developed within a model in which structures grow from small seeds, 
either Poisson distributed or set on a perturbed lattice. The extension of the PS work to more general and "standard" 
cosmological settings was due to Efstathiou, Fall & Hogan (1979), who limited their analysis to power-law spectra and an 
Einstein-de Sitter background. The first application of PS to a CDM spectrum was made by Schaeffer & Silk (1985). They 
"discovered' the PS MF (the procedure is described without any reference to the PS paper), complete with its unjustified 
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fudge factor 2, and criticized it in some interesting points; in particular, they tried to model, on purely geometrical grounds, 
objects which do not collapse spherically, like pancakes and filaments. 

Starting from 1988, many authors extended the PS approach in many directions, trying to understand why it appeared 
to work, in spite of its heuristic and not fully satisfactory derivation. The situation in those years is reviewed in Lucchin 
(1989). Just one year before, Kashlinsky (1987) tried to determine the 2-point correlation function of structures, collapsed 
according to the PS prescription. Martinez- Gonzalez & Sanz (1988a) performed similar calculations with a different 
method, while Martinez-Gonzalez & Sanz (1988b) attempted to determine the luminosity function of galaxies of different 
morphological types by means of an approach which was intermediate between the PS and the peak one. Lucchin & 
Matarrcsc (1988) formulated a PS MF for the case of non-Gaussian perturbations. Liljc (1992) and Lahav et al. (1991) 
extended the PS result to open Universes and to flat Universes with a cosmological constant. Zhan (1990) changed the 
PS procedure by taking into account the correction to the background density given by the initial density contrast; as a 
matter of fact, such a correction is negligible if the initial time is small. Schaeffer & Silk (1988a) used the PS procedure to 
justify the presence of some small-scale power in the HDM cosmology; in fact, the use of Gaussian smoothing causes some 
large-scale power to be spread toward small scales. Again Schaeffer & Silk (1988a, b), and Occhionero & Scaramella (1988) 
applied the PS formula to get many cosmological predictions on collapsed structures. BCEK solved the cloud-in-cloud 
problem using the "excursion set approach". Lacey & Cole (1994), introduced the merging histories concept, which gives 
an important piece of information in the formation history of dark matter objects. Del Popolo & Gambera (1998), ST, 
STl, ST2, Del Popolo (2002a), (2005), showed that the non-sphericity of collapse has important consequences on the MF. 

Press and Schechter were the first to performed N-body simulations to test the validity of their formula. They found 
some encouraging agreement, hut their simulations were limited to 1000 bodies, a very small number to reach any firm con- 
clusion. Efstathiou, Fall & Hogan (1979) performed similar simulations, with the same number of point masses, obtaining 
the same conclusions as PS. 

Later, Efstathiou et al. (1988) compared the results of larger (32^ P3M, scale-free power spectra) N-body simulations 
to the PS formula: their dynamical range in mass was large enough to test the knee of the MF. The surprising result was 
that the PS formula nicely fitted their abundances of simulated halos (as found by means of a percolation friend-of-friend 
algorithm). Further comparisons with N-body simulations were performed by Efstathiou & Rees (1988), Narayan & White 
(1988), Carlberg & Couchman (1989), Carlberg (1990), BCEK, Brainerd & Villumsen (1992), White, Efstathiou & Frenk 
(1993), Ma & Bertschinger (1994), Jain & Bertschinger (1994), Gelb & Bertschinger (1994), Katz, Quinn, Bertschinger & 
Gelb (1994), Lacey & Cole (1994), Efstathiou (1995), Klypin & Rhee (1994), Klypin et al. (1995), Bond & Myers (1996b), 
Governato et al. (1998), YNY. Most authors reported the PS formula to fit well their N-body results; nonetheless, all the 
authors agree in stating the validity of the PS formula to be only statistical, i.e. the existence of the single halos is not 
well predicted by the linear overdensity criterion of PS (see in particular BCEK). 

There are however some exceptions to this general agreement: Brainerd & Villumsen (1992) reported their MF, based 
on a GDM spectrum, to be very similar to a power-law with slope —2, different from the PS formula both at small and at 
large masses. Jain & Bertschinger (1994), Gelb & Bertschinger (1994) and Ma & Bertschinger (1994) noted that, to make 
the PS formula agree with their simulations (based on GDM or GHDM spectra), it is necessary to lower the value of the 
Sc parameter as redshift increases. The same thing was found by Klypin et al. (1995), but was interpreted as an artifact 
of their clump-finding algorithm. More recent simulations seem to confirm this trend (Governato et al. 1999; YNY)). 

Lacey & Cole (1994) extended the comparison to N-body simulations to the predictions for merging histories of dark- 
matter halos; they found again a good agreement between theory and simulations. This fact is noteworthy, as merging 
histories contain much more detailed information about hierarchical collapse. 

It is opportune to comment on two important technical points about such comparisons. First, the 6c parameters used 
by different authors as "best fit" values range from the 1.33 of Efstathiou & Rees (1988) to the 1.9 found (in a special case) 
by Lacey & Cole (1994). The precise value of the 6c parameter is influenced by the shape of the filter used to calculate 
the PS, Gaussian filters requiring lower 6c values. Recent simulations tend to give 6c — 1.5 (e.g. Klypin et al. 1995) or 
6c = 1.69 (e.g. Lacey & Cole 1994). If 6c changes with time, a value 1.5 could be good at high redshifts, lowering to 1.7 
at low redshifts. 

Second, the numerical MF deeply depends on the way halos are picked up from simulations. Typical algorithms, such 
as the friend-of-friend or DENMAX, are parametric, i.e. they contain free parameters. For instance, the frequently used 
friends-of-friends algorithm defines as structures those clumps whose particles are separated among them by distances 
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smaller than a percolation parameter b times the mean interparticle distance. A heuristic argument, based on spherical 
collapse, suggests a value of 0.2 for b (with this percolation parameter the mean density contrast of halos is about 180, 

which is the expected density contrast of a virialized top-hat perturbation) . Obviously, the use of different b parameters 
leads to different MFs. In practice, what is obtained in this case is not "the" MF, but the "fricnd-of-fricnd, 6=0.2" 
MF. Then, the numerical MF contains some hidden parameters, which, together with the dc parameter (and the mass 
associated to the filter in the Lacey & Cole (1994) paper), makes such comparisons more similar to parametric fits, rather 
than to comparisons of a theory to a numerical experiment. 

The present paper is organized as follows: section 2 contains some tools needed to formulate the MF theory: it is focused 
on simple models for structure formation (density perturbation growth, spherical collapse model, Zel'dovich approximation, 
ellipsoid collapse model, etc.). 

Section 3 contains a review of the theoretical MF problem and describes the theoretical works (such as the excursion 
set model) which have tried to extend the validity of the original Press & Schechter work, or have proposed alternative 
procedures. 

In the same section, it is reviewed the building of a MF theory, based on more realistic approximations for gravitational 
collapse than the spherical collapse. Finally, it is reviewed the work of Del Popolo & Gambera (1998), ST, STl, ST2, Del 
Popolo (2002a), (2005), who showed that the non-sphericity f collapse has important consequences on the MF. Section 4 
is devoted to prospects and conclusions. 



2. Theoretical bases of MF theory 

2.1. Background cosmology 

The simplest cosmological model that describes, in a sufficient coherent manner, the evolution of the universe, from 10~^s 
after the initial singularity to now, is the so called Standard Cosmological Model (or Hot Big Bang model). It is based 
upon the Priedmann-Robertson- Walker (FRW) metric, which is given by: 



ds' = c'de - a{tf 



1 — fcr^ 



(1) 



where c is the light velocity, a(t) a function of time, or a scale factor called "expansion parameter" , t is the time coordinate, 
r, 9 and (j) the comoving space coordinates. The evolution of the universe is described by the parameter a(t) and it is 
fundamentally connected to the value p of the average density. 

The equations that describe the dynamics of the universe are the Friedmann's equations (Friedman, 1924) that we are 
going to introduce in a while. These equations can be obtained starting from the equations of the gravitational field of 
Einstein (Einstein, A., 1915): 

„ 1 „ SttG 

Rik — -zQikii — J-Tik (2) 

2 c* 

where now, Rik is a symmetric tensor, also known as Ricci tensor, which describes the geometric properties of space-time, 
Qik is the metric tensor, R is the scalar curvature, Tik is the energy-momentum tensor. 

These equations connect the properties of space-time to the mass-energy. In other terms they describe how space-time is 
modeled by mass. Combining Einstein equations to the FRW metric leads to the dynamic equations for the expansion 
parameter, a(t). These last are the Friedmann equations: 

d{pa^) = -pd{a^) (3) 

1 .2 k 8ttG 
^"+^ = — ^ 

2- + ^ + 4 = -SttGp (5) 

where p is the pressure of the fluid of which the universe is constituted, k is the curvature parameter and a(t) is the 
scale factor connecting proper distances r to the comoving ones x through the relation r = a(t)x. Only two of the 
three Friedmann equations are independent, because the first connects density, p to the expansion parameter a(t). The 
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Fig. 1. Evolution of tlie scale factor. 

character of the solutions of these equations depends on the value of the curvature parameter, k, which is also determined 
by the initial conditions by means of Eq. 3. The solution to the equations now written shows that if p is larger than 
Pc = 1^ = 1.88 * 10~^^ g/cm^ (critical density, which can be obtained from Priedmann equations putting t = to, k = 0, 

and H = lOOkm/sMpc), space-time has a closed structure (k = 1) and equations shows that the system go through a 
singularity in a finite time. This means that the universe has an expansion phase until it reaches a maximum expansion 
after which it recoUapse. If p < pc, the expansion never stops and the universe is open fc = — 1 (the universe has a structure 
similar to that of an hyperboloid, in the two-dimensional case). If finally, p = pc the expansion is decelerated and has 
infinite duration in time, k = 0, and the universe is flat (as a plane in the two-dimensional case). The concept discussed 
can be expressed using the parameter = In this case, the condition = 1 corresponds to fc = 0, < 1 corresponds 
to fc = —1, and > 1 corresponds to fc = 1. Fig. 1 plots the evolution of the expansion parameter for k = —1, 0, 1. 
In a flat model, the FRW equation becomes: 

= -^-Gp, (6) 
whose solution is: 

a{t) = {t/tof'\ (7) 
In the case of open models with no cosmological constant: O < 1, A = 0, we can write: 

= |7rGp(l+(Oo-i-l)a), (8) 
and the a{t) evolution can be expressed through the following parametric representation: 

"^^^ = 2(l-°»o) ^'°'^^~^^ 

In the case of flat models with positive cosmological constant: ^l<l,A^0,f2 + A/3Hq = 1, we can write: 

^ = l^Gp{l+{n,'-l)a'), (10) 



a(f) = (l7„i-l)-^/%inhV3 . (11) 



d2 




The value of Cl can be calculated in several ways. The most common methods are the dynamical methods, in which 
the effects of gravity are used, and kinematics methods sensible to the evolution of the scale factor and to the space-time 
geometry. 



2.2. Perturbations evolution. 

Density perturbations in the components of the universe evolve with time. In order to get the evolution equations for S 
in Newtonian regime, it is possible to use several models. In our model, we assume that gravitation dominates on the 
other interactions and that particles (representing galaxies, etc.) move coUisionless in the potential of a smooth density 
function (Peebles 1980). 

The distribution function of particles for position and momentum is given by: 



dN = /(x, p, t)(fx(fp 



(12) 
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and density: 

p(x, t) = ma-^ J p, t) ^pb[l + (5(x, t)] (13) 

where m is the mass of a particle and pb the background density. Applying Liouville theorem to the probability density 
on a limited region of phase-space of the system we have that f verifies the equation: 

§7 + -^V/-mV'^|^=0 (14) 
ot ma-' op 

The distribution function f that appears in the previous equations cannot be obtained from observations. It is possible 
to measure moments of f (density, average velocity, etc.). We want now to obtain the evolution equations for S. For this 
reason, we start integrating Eq. (|14|l on p and after using Eq. H13|l . we get: 

or -I p 

a'pt^ + ^\7 Pfd'p^O (15) 



(16) 



dt 

If we define velocity as: 

J ma '' ^ 

Jfd^p 

and introduce it in Eq. I|15f) we get: 

Pb^ + -V (pv) = (17) 
ot a 

We can now multiply Eq. 114() for p and integrate it on the momentum: 
d f If 

Pafd^p + ^dp / PcPpfd^p + a^p(x, i)0,„ = (18) 



dt 

this last in Eq. H15() leaves us with: 

^ + 2^- . - V [(1 + ^) V + ^d^d, j p^p.H'p (19) 
and finally using 

^ J fPaPpd^P 
ma'' j jd'^p 

the equation for the evolution of overdensity becomes: 

£ + 4i = i ^ + ^) ^ + ^^"^'^ + '^^ < ^^^^ 

(Peebles 1980). The term < VaVp > is the tensor of anisotropy of peculiar velocity. This is present in the gradient and 
then it behaves as a pressure force. If we consider an isolated and spherical perturbation, it is possible to assume that 
initial asymmetries does not grow up and so we can suppose, in this hypothesis that < VaVfj >= 0. 

When the density contrast is much smaller than one, ^ ^ 1, and peculiar velocities, v, are small enough to satisfy 
(vt/d)'^ ^ S, where t is the cosmological time and d is the coherence length of the matter field, one can obtain a linear 
theory for perturbations evolution as follows. 

In this case we have: 

This equation in an Einstein-de Sitter universe (fi = 1, A = 0) has the solutions: 

S+ = A+(x)ti 5-(x,i) = A-{x)t-'^ (23) 

The perturbation is then done of two parts: a growing one (which shall be denoted with b{t)), becoming more and more 
important with time, and a decaying one becoming negligible with increasing time, with respect to the growing one. 
The solutions for the growing modes, relative to the three background models, are reported: 
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= 1: 

b{t) = a{t). (24) 
f2 < 1, A = 0: it is useful to use the time variable 



r = (1 - mr'^' = V («("o ' - 1))"' + 1- (25) 

Then: 

Kr) = ^,J, f 1 + 3(r^ - 1) fl + I In f ) ) . (26) 



2(^^0-1 - 1) V ' V 2 

Note that this b{t) function saturates to the value 5/2(^1,^^ — 1) at large times, 
n < 1, A ^ 0, flat: it is useful to use the time variable 

h = coth(3tVA73/2). (27) 

Then, 

6(t) = h {x'^ix^ - iy/^)-^dx. (28) 

Jh 

Growing modes are normalized so as to give b{t) ~ a{t) at early times, and a{to) = 1. 

In the MF theory, collapse time estimates are often based on an extrapolations of the linear regime to density contrasts 
of order one. It is then convenient to define the quantity: 

Si = 5{U)/h{U). (29) 

This is the initial density contrast linearly extrapolated to the time at which h{t) = 1, which, in an Einstein-de Sitter 
background, is the present time; it will be used in the next sections. 

Before concluding this section, we want to find an expression for the velocity field in the linear regime. Using the 
equation of motion p = ma^x, ^ = — m V (t> and the proper velocity of a particle, v = ax, verify the equation: 

^ + = = Gp,a [ d^xSi.,t)^^ (30) 

at a a J |x — x'| 

Supposing that v is a similar solution for the density, v = V+(x, t)t^, we get: 

= I d'-'Pr^ (31) 

4tt dx" y |x' - x| ^ ^ 

(Peebles 1980). This solution is valid just as that for 6 in the linear regime. At time t = to this regime is valid on scales 
larger than 8h~^Mpc. 



2.3. The spectrum of density perturbation. 

In order to study the distribution of matter density in the universe it is generally assumed that this distribution is given 
by the superposition of plane waves independently evolving, at least until they are in the linear regime (this means till 
the overdensity 5 = << 1). Let we divide universe in cells of volume Vu and let we impose periodic conditions on the 
SLufaces. If we indicate with p the average density in the volume and with p(y) the density in r, it is possible to define the 

density contrast as: 

5{v) = (32) 
This quantity can be developed in Fourier series: 

5{r) = (5kea;p(ikr) = 6\^exp{—ik.Y) (33) 

k k 
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( Kolb e Turner 1990), where kx = ^^"^ (and similar conditions for the other components) and for the periodicity condition 
S{x, y, L) = d{x, y, 0) (and similar conditions for the other components). Fourier coefficients (5k are complex quantities given 
by: 

Sk = ^ [ S{r)exp{-ikr)dr (34) 
III Jv„ 

For mass conservation in Vu we have also (5k=o = while for reality of S{r), (5{^ — (5_k. If we consider n volumes, Vu, we 
have the problem of determining the distribution of Fourier coefficients 6k and that of \5\. We know that the coefficients 
are complex quantities and then 6^ = \6k \ exp{i9k). If we suppose that phases are random, in the limit y„ — > oo it is 
possible to show that we get \d\^ = J^w I'^kl^- The Central limit theorem leads us to conclude that the distribution for 5 
is Gaussian: 

P{S) cx exp{—) (35) 
(Efstathiou 1990). The quantity a that is present in Eq. H35|) is the variance of the density field and is defined as: 

k " k 

This quantity characterizes the amplitude of the inhomogeneity of the density field. If Vu ^ oo, we obtain the more usual 
relation: 

^ ' P{k)(f-k = ^ f P{k)k^dk (37) 



{2TTf J ' ' 



The term P{k) =< \d\ > is called "Spectr um of perturbations" . It is function only of k because the ensemble average 
in an isotropic universe depends only on r. A choice often made for the primordial spectrum is P{k) =^ Ak^^ which in the 
case n = 1 gives the scale invariant spectrum of Harrison- Zeldovich. An important quantity connected with the spectrum 
is the two-points correlation function ^(r,t). It can be defined as the joint probability of finding an overdensity S in two 
distinct points of space: 

e(r,i) (5(r,t)5(r + x,i) > (38) 

(Peebles 1980), where averages are averages on an ensemble obtained from several realizations of universe. Correlation 
function can be expressed as the joint probability of finding a galaxy in a volume SVi and another in a volume SV2 
separated by a distance ri2: 

S^P = n^[l + ^{ri2)]SViSV2 (39) 

where ny is the average number of galaxies per unit volume. The concept of correlation function, given in this terms, can 
be enlarged to the case of three or more points. 

Correlation functions have a fundamental role in the study of clustering of matter. If we want to use this function for 
a complete description of clustering, one needs to know the correlation functions of order larger than two (Fry 1982). 
By means of correlation functions it is possible to study the evolution of clustering. The correlation functions are, in 
fact, connected one another by means of an infinite system of equations obtained from moments of Boltzmann equation 
which constitutes the BBGKY (Bogolyubov-Green-Kirkwood-Yvon) hierarchy (Davis e Peebles 1978). This hierarchy can 
be transformed into a closed system of equation using closure conditions. Solving the system one gets information on 
correlation functions. 

In order to show the relation between perturbation spectrum and two-points correlation function, we introduce in Eq. 
H38|l . Eq. recalling that 5^ = <5(— k) and taking the limit Vu — > 00, the average in the Eq. H38(l can be expressed in 

terms of the integral: 

m = I \S{k)fexp{-tkr)d'k (40) 

This result shows that the two-point correlation function is the Fourier transform of the spectrum. In an isotropic universe, 
it is |r| = r and then |k| = k and the spectrum can be obtained from an integral on |k| = k. Then correlation function 
may be written as: 

ar) = mkf-^dk (41) 
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During the evolution of the universe and after perturbations enter the horizon, the spectrum is subject to modulations 
because of physical processes characteristic of the model itself (Silk damping (Silk 1968) for acollisional components, free 
streaming for collisional particles, etc.). These effects are taken into account by means of the transfer function T{k;t) 
which connects the primordial spectrum P(fc; tp) at time tp to the final time tf. 

2 

T^{k-tf)P{k-tp) (42) 

where b(t) is the law of grow of perturbations, in the Hnear regime. In the case of CDM models the transfer function is: 

-1 

T{k) = + [ak + [bkf-'^ + [ckfY^ " (43) 

(Bond e Efstathiou 1984), where a = &A{Vlh'^)-^Mpc; b = 3.0(f]/i2)"iAfpc; c = 1.7{nh'^)-'^Mpc; v = 1.13. It is interesting 
to note that Eq. (|35|l is valid only if tr << 1, since \b\ < 1. This implies than non- linear perturbations, ct >> 1, must be 
non-Gaussian. In fact when the amplitudes of fluctuations grow up, at a certain point modes are no longer independent 
and start to couple giving rise to non-linear effects that change the spectrum and correlation function (Juskiweicz et al 
1984 ). There are also some theories (e.g., cosmic strings (Kibble & Turok 1986)) in which even in the linear regime 
perturbations are not Gaussian. 



P{k;tf) = 



Htf) 



2.4. Spherical Collapse 

Linear evolution is valid only if ^ << 1 or similarly, if the mass variance, a, is much less than unity. When this condition 
is no longer verified (e.g., if we consider scales smaller than 8/i~^ Mpc), it is necessary to develope a non-linear theory. 
In regions smaller than Sh~^ Mpc galaxies are not a Poisson distribution but they tend to cluster. If one wants to study 
the properties of galactic structures or clusters of galaxies, it is necessary to introduce a non-linear theory of clustering. A 
theory of this last item is too complicated to be developed in a purely theoretical fashion. The problem can be faced, by 
using N-Body simulations of the interesting system, assuming certain approximations that simplifies it. Spherical Collapse 
Model, Zel'dovich approximation (Zel'dovich 1970). This last gives a solution to the problem of the grow of perturbations 
in an universe with p = not only in the linear regime but even in the mildly non-linear regime. 

Spherical symmetry is one of the few cases in which gravitational collapse can be solved exactly (Gunn & Gott 1972; 
Peebles 1980). In fact, as a consequence of Birkhoff's theorem, a spherical perturbation evolves as a FRW Universe with 
density equal to the mean density inside the perturbation. 

The simplest spherical perturbation is the top-hat one, i.e. a constant overdensity 5 inside a sphere of radius i?; to 
avoid a feedback reaction on the background model, the overdensity has to be surrounded by a spherical underdense shell, 
such to make the total perturbation vanish. The evolution of the radius of the perturbation is then given by a Friedmann 
equation. 

The evolution of a spherical perturbation depends only on its initial overdensity. In an Einstein-de Sitter background, 
any spherical overdensity reaches a singularity (collapse) at a final time: 

= y (^<5(t.)) "'u- (44) 
By that time its linear density contrast reaches the value: 

5i{tc) = 5, = l^(^-^^ ' - 1.69. (45) 

In an open Universe not any overdensity is going to collapse: the initial density contrast has to be such that the total 
density inside the perturbation overcomes the critical density. So, a perturbation ought to satisfy 5i > 1.69 • '2,{riQ^ — l)/5 
to be able collapse. 

Of course, collapse to a singularity is not what really happens in reality. It is typically supposed that the structure 
reaches virial equilibrium at that time. In this case, arguments based on the virial theorem and on energy conservation 
show that the structure reaches a radius equal to half its maximum expansion radius, and a density contrast of about 
178. In the subsequent evolution the radius and the physical density of the virialized structure remains constant, and its 
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density contrast grows with time, as the background density decays. Similarly, structures which collapse before are denser 
than the ones which collapse later. 

Spherical collapse is not a realistic description of the formation of real structures; however, it has been shown (see 
Bernardeau 1994 for a rigorous proof or Valagcas 2002a, b) that high peaks (> 2a) follow spherical collapse, at least in the 
first phases of their evolution. However, a small systematic departure from spherical collapse can change the statistical 
properties of collapse times. 

Spherical collapse can describe the evolution of underdensities. A spherical underdensity is not able to collapse (unless 

the Universe is closed!), but behaves as an open Universe, always expanding unless its borders collide with neighboring 
regions. At variance with overdensities, underdensities tend to be more spherical as they evolve, so that the spherical 
model provides a very good approximation for their evolution. 



2.5. Zel'dovich approximation. 

Most MF theories proposed in the literature are based at best on spherical collapse, which neglects tides which are the 
relevant dynamical interaction. Spherical top-hat collapse is a truly local dynamical approximation: the fate of a spherical 
perturbation is determined just by its initial overdensity. In other words, the dynamical role of the whole Universe outside 
the perturbation (according to Birkhoff 's theorem) is assumed to be negligible. It is possible to construct a mixed Eulerian- 
Lagrangian system from the evolution equations of fluid elements by decomposing the tensor of (Eulerian) space derivatives 
of the peculiar velocity u into an expansion scalar 0, a shear tensor aab and a vorticity tensor ujah. In this way, the following 
evolution equation for the density contrast can be obtained (see, e.g., Ellis 1971; here the growing mode b{t) is used as 
time variable): 

Here = <Tab<^ab/'^ and = uiab^^ab/^ (note that in this context is not the mass variance). According to linear 
theory, any vortical mode is severely damped in the early gravitational evolution, and this remains true during the mildly 
non- linear regime, up to "orbit crossing OC ^ (at OC vorticity couples with the growing mode; see Buchert 1992). Then 
it is reasonable to assume vanishing vorticity in the present framework. On the other hand, the shear does not vanish in 
general: it provides the link between the mass element and the rest of the Universe. 
The evolution equation for the shear reads as follows: 

+ ^'^cTab + (Tacf^cb + ^nGp^aab - ^cr'^Sab = -^irGp^Eab- (47) 

The tensor Eab, represents the tidal interactions between the mass element and the rest of the Universe. Then, tides are 
the relevant dynamical interaction neglected by spherical collapse. 

It is useful to find which is the simplest way to introduce tides in the evolution of a mass element. The simplest, 
realistic approximation of gravitational evolution in the (mildly) non-linear regime is the Zel'dovich approximation, that 
I am going to summarize. 

In this approximation, one supposes to have particles with initial position given in Lagrangian coordinates q. The 
positions of particles, at a given time t, the Lagrangian-to-Eulerian mapping is written as follows: 

x = q + 6(t)p(q) (48) 

where x indicates the Eulerian coordinates, p(q) describes the initial density fluctuations and b{t) describes their grow in 
the linear phase and it satisfies the equation: 

— +2a-^ — — = AnCpb (49) 
dt^ dtdt'^ ^ ' 

The equation of motion of particles, according to the quoted approximation, is given by: 

V = dq + 6p(q) (50) 

^ In the dynamical evolution of cold matter, it can happen that two mass elements get to the same point. This event is called 
"orbit crossing" 
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The peculiar velocity of particles is given by: 

while the density of the perturbed system is given by: 





dqj 




p 


dxk 


= P 



(52) 



Developing the Jacobian present in Eq. 152|l at first order in b{t)p{q), one obtains: 



y «-6(t)VqP(q) (53) 

This equation can be re-written, separating the space and time dependence, as in the equation for u, and writing: 

b{t) = ti p(q) = ^i-^Akea;p(?kq) (54) 

in the form: 

^ = A]^t^ expiikq) (55) 
^ k 

(Efstathiou 1990), that leads us back to the linear theory. In other words, Ze'ldovich approximation is able to reproduce 
the linear theory, and is also able to give a good approximation in regions with ^ >> 1. Using the expression for p{q), 
the Jacobian in Eq. 1)52(1 is a real matrix and symmetric that can be diagonalized. With this p(q) the perturbed density 
can be written as: 

" (1 - 6(t)Ai(q))(l - b{t)X,{q)){l - bit)Xs{q)) ^^^^ 

where Ai, A2, A3 are the three eigenvalues of the Jacobian, describing the expansion and contraction of mass along the 
principal axes. From the structure of the last equation, we notice that in regions of high density Eq. H5t)|l becomes infinite 
and the structure of collapse in a pancake, in a filamentary structure or in a node, according to values of eigenvalues. 
Some N-body simulations (Efstathiou e Silk, 1983) tried to verify the prediction of Ze'ldovich approximation, using initial 
conditions generated using a spectrum with a cut-off at low frequencies. The results showed a good agreement between 
theory and simulations, for the initial phases of the evolution {a{t) ~ 3.6). Going on, the approximation is no more valid 
starting from the time of shell-crossing. After shell-crossing, particles does not oscillate any longer around the structure 
but they pass through it making it vanish. This problem has been partly solved supposing that particles, before reaching 
the singularity they sticks the one on the other, due to a dissipative term that simulates gravity and then collects on the 
forming structure. This model is known as "adesion- model" (Gurbatov et al. 1985). 

Summarizing, Zel'dovich approximation gives a description of the transition between linear and non-linear phase. It is 
expecially used to get the initial conditions for N-body simulations. The main problem with Zel'dovich and perturbative 
Lagrangian approximations is that they break down after OC. Many authors have then tried to develop approximations 
which avoid OC or make particles oscillate around pancakes (see Sahni & Coles 1996 for a complete review). It is interesting 
to note that the Zel'dovich approximation is the first term of a perturbative series; the perturbed quantity is not the density, 
as in Eulerian perturbation theory, but the displacement of the particles from the initial position. 

2.5.1. Collapse time in Zel'dovich approximation 

In the Ze'ldovich approximation, the density contrast b, the expansion and the shear Oab evolve as follows: 

J(q, i) = (1 - fe(i)Ai)(l - b{i)\^){\ - b{t)X3) (57) 
Al A2 A3 



1 - b{t)Xi 1 - b{t)X2 1 - b{t)X3 
Al 9 X2 6 



a^,{q,t) ^ diag _ ^^^^^^ ^, ^ _ ^^^^^^ ^ _ ^^^^^^ 
5(q, t) = ((1 - 6(0Ai)(l - b{t)X2){l - b{t)X3)y 
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When b{t) — 1/Ai, caustic formation takes place: the Jacobian determinant vanishes, and all the other quantities go to 
infinity. It is then quite reasonable, from the point of view of the mass element, to define such instant as the collapse time. 
Hereafter collapse will always be defined as the OC event. Note also that, after OC, the Zel'dovich evolution along the 
second and third axes is not really meaningful, as Zel'dovich does not work after OC. It is then possible to give collapse 
time estimates for any mass element. Initial conditions are "locally" given by the three A eigenvalues, but the evolution is 
not physically local, as initial conditions contain non-local information about tides 

It is useful to define the following variables (M98): 



a; = Ai - A2 (58) 
y = A2 - A3, 

and to use the growing mode b(t) as time variable. Moreover, it is possible to consider regions with linear initial density 
contrast 5/ = 1 or —1, as all the other cases can be obtained by a simple rescaling of b. The Zel'dovich collapse time is 
finally given by: 



el' dovich 



5i + 2x + y 



(59) 



Fig. 2 shows the collapse time curves bc{x,y) for Si — 1 and —1. A problem is soon apparent: in the spherical case, 
when X = y = Q, the collapse time is 3 /Si, instead of the exact 1.69/(5; value. This discrepancy is due to the fact that the 
Zel'dovich approximation is an exact solution (before OC of course) in one dimension (see, e.g., Shandarin & Zel'dovich 
1989), and is then able to describe the collapse of pancake- like structures, while it severely underestimates the collapse 
rate in spherical symmetry. 

A way to overcome this problem is to try some simple ansatze for the "true" shape of the collapse time curve. In 
practice, a truly realistic collapse time curve will depend not only on the A eigenvalues, but also on the values of the 
density (or potential) field in all the points of (Lagrangian) space. 

A first way to change the Zel'dovich prediction is to force it not to assume values larger than the spherical one (M98): 

= min{6f 1.69/(5,}. (60) 

This be curve is shown in Fig. 3a for Si = 1; it has a plateau, of height 1.69, at small x and y values. On the other 
hand, it is unlikely that all quasi-spherical collapses have exactly the same be value; a systematic trend of lower b^. values 
at nonvanishing x and y values is more realistic. Such a be curve can be modeled as the intersection of the Zel'dovich 
prediction with a slightly inclined plane which reaches 1.69 at a; = y = (M98): 

- min{6f 1.69/(5, - e{2x + y)}. (61) 

This curve, with e — 0.2, is shown in Fig. 3b. Monaco (1995) contains a more complete discussion of such ansatze. A further 
possibility, examined in Monaco (1995), is to insert "by hand" in Eq. (|46|l the shear as given by Zel'dovich fEa. I57|) . and 
then solve the equation numerically; the resulting collapse time is very similar to that of Fig. 3a, but the transition from 
spherical to Zel'dovich regime is smooth. 

An interesting conclusion, is that the reasonable systematic trend shown in Fig. 3b does influence the large-mass part 
of the mass function, moving it toward large masses, even though spherical collapse is recovered for very large overdensities 
(which are characterized by small x and y values). In other words, the fact that large, rare fluctuations asymptotically 
follow spherical collapse does not guarantee that the "spherical" PS MF (with Sc = 1.69) is recovered at large masses. 



2.6. Ellipsoidal collapse 

The convenience in using the homogeneous ellipsoid collapse model resides in the fact that it can easily be solved by 
means of a numerical integration of a system of three second-order ordinary differential equations. One of the advantages of 
spherical symmetry is that, because of Birkhoff 's theorem, it is possible to introduce in a background metric a perturbation 
without influence the rest of the Universe, provided any positive perturbation is compensated for by an (outer) negative 
one, such to make the total mass perturbation vanish. This is necessary to ensure the self-consistency of the problem: 
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Fig. 2. Collapse times with Zel'dovich approximation. Figures taken by M98. 




Fig. 3. Collapse times with the two ansatze. Figures taken by M98 



the background has to evolve as if it were unperturbed. This reasoning is not valid any more when introducing a triaxial 
perturbation in an luipcrturbcd background: this is going to influence the background, through non-linear feedback effects. 
To use ellipsoidal collapse in a cosmological context, the correct strategy is not to try to insert an ellipsoid in a uniform 
background, but to extract an ellipsoid from a perturbed FRW Universe. 

In the following, we shall describe the ellipsoidal model and how to find solutions to the equations, both numerically 
and analytically. 

The dynamical variables of ellipsoidal collapse are the three axes ai{t) of the ellipsoid; they are normalized as the scale 
factor: ai{t) = a{t) if the ellipsoid is a sphere with null density contrast. Their evolution equations are: 



^ - (2a(l + in,^ - l)a))-'^ + (2a2(l + [n,' l)a))-'a, 



3 + 3+ 2^ + A„, 



= 



(62) 



in the open case (the Einstein-de Sitter case can be obtained by setting fio = 1)) while in the flat case with cosmological 
constant they are: 



(fai _ 1 - 2(^0^ - l)a^ dui 
~da? 2a(l + {^q^ - l)a^) da 



+ (2a2(l + (Oq 1 - l)a))-^ai 



0. 



(63) 
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Note that this time the scale factor a{t) has been used as time variable. The density contrast 6 is: 

5 = - 1, (64) 

while the quantities b'^ and A^^ are defined as: 
2 

K = g [ajajOfc^Dla- , a^, 4) - 1] i^j^k (65) 
(where the is the Carlson's elliptical integral 

Rn{x,y,z) = -j^ (, + ,)i/2(,+^)i/2(, + ,)3/2 > (66) 
which can be calculated by means of the routine given by Press & Teukolsky (1990)) and 

A;, = -^(^-aoA.). (67) 

Initial conditions can be set by imposing that the evolve according to Zel'dovich approximation at early times: 

Ui ~ a(l — aXi) (68) 

^~i(a,(a)-a2A,). (69) 
da a 

These three coupled second-order ordinary differential equations, Eqs. (|62|) or Ht)3|l . can be integrated by means of 
standard routines, as the Runge-Kutta one given in Press et al. (1992). The numerical integration has to be pushed to the 
singularity, when at least one axis vanishes (and the density diverges). To do this, it is useful to use logarithmic variables, 
to have more controlled variations from quasi-homogeneity to collapse. Moreover, the integration can be divided into two 
parts: the first is stopped at decoupling, defined as the instant at which the density starts to increase, while in the second 
part the density is used as time variable, and the integration is pushed up to large density values (Monaco 1997a). The 
overall precision of the numerical integration is better than 1% for the spherical collapse, but becomes about 8% for 
pancake-like collapses. It is to be noted that, in any case, the first collapsing axis is that corresponding to Ai, the largest 
A eigenvalue. 

Fig. 4 shows the collapse "times" be of initially overdense {5i=l) and underdense ((5;=— 1) ellipsoids in an Einstein- 
de Sitter Universe (in this case be = ftc). Spherical collapse is obviously recovered ai x ^ y ~ 0, while quasi-spherical 
collapses reasonably show a systematic departure from spherical collapse, as in Fig. 3a. The large-shear behavior is similar 
but not identical to that predicted by Ze'ldovich approximation (ZEL); at variance with what happens with quasi-spherical 
collapses, in this range ZEL tends to underestimate the collapse time. Fig. 5 shows the b^ curve for ellipsoids in an open 
Universe; this time 6i=3 has been chosen, to allow all the ellipsoids to collapse. As expected, this curve is nearly identical 
to the one shown in Fig. 3a, apart from an obvious rescaling. Notably, numerical calculations of collapsing ellipsoids in 
open Universes are affected by larger errors than that quoted above if the ellipsoid takes a long time to collapse. 

An analytic solution to the collapse model equations can be obtained as follows. For the unisolated ellipsoid, the 
evolution equations reduces to three equations for the three semiaxes of the ellipsoid and are given by (Watanabe 1993; 
van de Weygaert 1996): 



^ = -27rG \ po(ai - 7&i) + 



2^G7 (-6i) (pe - Ph) ai (70) 



where: 
3 Q 



b=(-/3,/3-l,l) (71) 

ZTT 
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Fig. 4. Collapse times with Ellipsoidal model. Figures taken by M98. 




Fig. 5. Collapse times with Ellipsoidal model, open Universe. Figures taken by M98. 



and where pb is the density of the background universe, and pc the density within the ellipsoid. The coefficients ai are 
given by: 

dX 



tti = aia2a3 / (72) 

Jo {af + X)[{al + X){al + X)ial + X)]^ 

Assuming that the external structures, giving rise to the tidal field, are at a large distance from the ellipsoid (see 
Eisenstein & Loeb 1995), the amplitude of the external quadrupole force is assumed to increase with the linear growth 
rate (Ryden 1988; Watanabe 1993; Eisenstein & Loeb 1995), D{t) (this last quantity is given in Peebles 1980): 
D{t) 



Q{t) = Qo- 



Dn 



(73) 



Here the subscript "0" means that the corresponding quantity is calculated at the present epoch and D{t) — Rh{t), for an 

flb(to) 



Einstein-de Sitter (hereafter EdS) universe. In an Einstein-de Sitter universe, Eq. CZai, Q(0 = Qo o^Vi reduces to 



Q{to) = Qo time to. 

In order to have an estimate of the value of Qo, for a cluster interacting with a neighboring one, one can use the simple 
model in Watanabe (1993), considering a cluster which has a neighboring cluster with a mean density contrast < 5 >~ 3, 
a comoving separation (0,0, X3), and a comoving size Axa — X3/3, the Q33 quadrupole component is given by: 



Q 



33 



8 _ / Ax3 



0.3 



(74) 



The previous estimate corresponds to a cluster interacting with a neighbor having a mass excess comparable to that of 
the Virgo cluster, and a separation three times its size. 
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Unisolated ellipsoid 

Numerical solution 

Analytical solution 

Axial ratio: 1:1.25:1.5 

Qo=0.1; b. = {-0.5,-0.5,l) 




(a) 




Unisolated ellipsoid 

Numerical solution 

Analytical solution 

Axial ratio: 1:1.25:1.5 

Qo=0.3; b.=(-0.5,-0.5,l) 



(b) 




Unisolated ellipsoid 

Numerical solution 

Analytical solution 

Axial ratio: 1:1.25:1.5 

Q„=0.4; b, = ( 0.5, 0.5,1) 



500 1000 1500 

(c) 

Fig. 6. Evolution of unisolated homogeneous ellipsoidal perturbations in an EdS universe with Hq — 50km/s/Mpc, 
Pc/ph = 1.003, axial ratio is 1 : 1.25 : 1.5, k = (-0.5, -0.5, 1), and Qo = 0.1 (a), Qo = 0.2 (b), Qo = 0.4 (c). Figures taken 
from Del Popolo (2002b). 



Another way of estimating Qq is by using the anisotropy of the velocity field in the LSC from data of Lilje, Yahil & 
Jones (1986). If we indicate with Qvo the component of the largest absolute value of the anisotropic velocity, one gets: 
Qol7g-6 = ^Q^o (Watanabe 1993). Since Lilje, Yahil & Jones (1986) deduced a value of Qvo ~ 0.1 — 0.2 at the distance 
of the Local Group from Virgo, we have that Qq^Iq-^ ^ 0.4 — 0.8. 

It is possible to obtain an analytical solution of Eq. H7U|) . describing the evolution of an unisolated ellipsoid as shown 
in Del Popolo (2002b) In this case the solution can be written in the form: 

^=i?b-^ai(i?b-i?e)-dxi?l'^) fl-^) (75) 



^ = i?b - (i?b - Ro) (76) 

0-2 [ti ) ^ 
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(a) 



Numerical solution 

Analytical solution 

Axial ratio: 1:1.25:1.5 

Top line: %-0.2; bj-(-0.5,-Q.5, l) 
Medium line: Qq-O,!; b -(-0.5,-0.5, 1) 
Bottom line: Isolated spheroid 




(b) 



Numerical solution 

Analytical solution 

Axial ratio: 1:2:3 

Top line: Q„-0.2: b,-(-0.5,-0.5, 1) 
Medium line: Q„-0-l: b -(-0.5,-0-5, 1) 
Bottom line: Isolated spheroid 




Fig. 7. The evolution of the density contrast. The axial ratio of the ellipsoid is 1 : 1.25 : 1.5 (a) and 1:2:3 (b), the lines 
from bottom to top represent the case of an isolated ellipsoid (6i = (0, 0, 0)), Qo — 0.1, fci = (—0.5, —0.5, 1), and Qq — 0.2, 
bi ~ (—0.5, —0.5, 1), respectively. Figure taken from Del Popolo (2002b). 



7TT — — 7:"^3 (^b — Rc 

a^(ti) 2 



(77) 



where Rh is the scale factor of the background universe, and Rc that of the ellipsoid. 

For ellipsoids having the initial axial ratio 1 : a2 : as, with ai > 1.25 and a2 > 1.5, we now have that ci = 1.23, 
d — Q-K 10^'' and: 

0.15 



0.0672 ( ?i 1 U 



ai — ai 



^ — „0.07„-0. 06^-0. 01 



1,0.6 



02 + 0.031 



bi 



0.5 



as = 1.002aVa2o 



0.1„-0.035„-0.065 



lO.95 
"3 



0.09i,0.95\ 



-'SO 



as - 0.063ay^6^ 



(78) 

(79) 
(80) 



where bi was defined in Eq. H71|l . 

In the case of prolate spheroids, with axial ratio 1 : 1 : as and 1 < as < 5, a better approximation to the ai is: 



ai=ai + 0.037 ( — 
,aio 



, 0.35 /, X 0.15 
«30\ fh\ ^0.6 

b-?, J 



a2 



oio 
aso 



aso/ \bi 



0.5 



lO.95 
"3 



53 ^ 1.002 ( ^ ) (as - 0M3al °X 



0.065 



,0.09i,0.95\ 



(81) 
(82) 
(83) 



In the case of an unisolated ellipsoid the length of the uncoUapsed axes at collapse can be obtained by means of Eq. 
((7^ -((77 |) : 

as(ic) _ as(ti) ^b — §"3 (^b — -Rc) 



a2(ic) a2(ii)i?b- |52(i?b-i?c) 

The evolution of the density contrast can be calculated using the usual definition: 

g _ Pc - Ph _ PcO OlO 020 «30 / Rh \ ^ _^ 
Ph 



Pho ai a2 as \i?bo 



(84) 



(85) 
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12 3 4 



Fig. 8. The density contrast at turnaround for a prolate spheroid for several values of the longest axis, as, (the other 
two axes have fixed value ai : 02 = 1 : 1). The solid lines, from top to bottom, represent numerical results for the density 
contrast for an isolated spheroid (6i = (0,0,0)), and for unisolated spheroids with Qo = 0.1, bi = (—0.5,-0.5,1) and 
Qo = 0.2, (—0.5, —0.5, 1), respectively. The dashed lines represent the approximate solution (Eq. (37)). The upper dotted 
line represents the value of the density contrast at turnaround for a spherical perturbation. Figure taken from Del Popolo 
(2002b). 




2 4 6 



Fig. 9. Density contrast at virialization. The solid line refers to an isolated prolate spheroid. Figure taken from Del Popolo 
(2002b). 



If x{t) = XoX{t), y{t) = yoYit) and z{t) = ZoZ{t), are the principal axes (xq, yo and Zo are the initial values of the 
axes), the overdensity of the ellipsoid is the same used previously, the initial conditions are X = Y = Z = Rh = Rg = lat 
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Fig. 10. Turnaround epoch for a prolate spheroid. The sohd, short-dashed and fong-dashed hues, represent respectively 
the time of turnaround for an isolated spheroid and for spheroids having Qq — 0.1, bi — (—0.5,-0.5,1) and Qq = 
0.2, (-0.5,-0.5,1). The upper dotted line represents the value of the density contrast at turnaround for a spherical 
perturbation. Figure taken from Del Popolo (2002b). 

t = to and as before the initial velocity is equal to the Hubble velocity at to (representing the initial time) . The parametric 
equations satisfied by Rc(t) are: 

i?e = ^ (1 - cosm , ^ - ^ - sin(^)) (86) 

while, since our background is an EdS universe, Rh{t) ex t'^^^. 
It is easy to find that the density contrast is given by: 



XYZ 



-dll- 



2 / 6f- 



1 (87) 



where / = fii''^) = sin(z?)) and /2(i9) = 1 — cos('i9). The density contrast at turn-around is obtained by 

calculating (5i(t?ta), where the parameter t9ta at the turnaround epoch is given solving the equation: 



2 _ 4fi?r + if^/f - 1 



(88) 



35i dfRl^^ - 1 

Eq. H87() yields the familiar value S = (37r/4)^ in the spherical case, (rf ^ 0, ai = ai = a2 = a2 — cts = ois = 2/3). In 
general, in order to obtain the density contrast at turnaround, one has first to solve Eq. H88|l for for an arbitrary axial 
ratio and substitute the value in Eq. I|87|l . The time of turn- around can be calculated by: 

(i9 - sin(z?)) = ^ (i9 - sin(i9)) (89) 



4(53/2 ^ ' 2tt 
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where tg is the free-fall time: 

= (90) 
Similarly it is possible to calculate the density contrast at collapse and collapse time. 

Figs. 6-10 plots the comparison of the analytical model with numerical simulations for several interesting quantities. 

Fig. 6 shows the evolution of unisolated homogeneous ellipsoidal perturbations in an EdS universe with _ffo = 
50km/s/Mpc, pe/pb = 1-003, axial ratio is 1 : 1.25 : 1.5, h = (-0.5, -0.5, 1), and Qo = 0.1 (a), Qo = 0.2 (b), Qo = 0.4 (c). 

Fig. 7, plots the evolution of the density contrast. The axial ratio of the ellipsoid is 1 : 1.25 : 1.5 (a) and 1:2:3 (b), 
the lines from bottom to top represent the case of an isolated ellipsoid (6i = (0,0,0)), Qo = 0.1, bi = (—0.5, —0.5, 1), and 
Qo = 0.2, 6i = (—0.5, —0.5, 1), respectively. 

Fig. 8, shows the density contrast at turnaround for a prolate spheroid for several values of the longest axis, 03, (the 
other two axes have fixed value ai : 02 = 1 : 1). The solid lines, from top to bottom, represent numerical results for the 
density contrast for an isolated spheroid {bi = (0,0,0)), and for wmsotoerf spheroids with Qo = 0.1, bi = (-0.5,-0.5,1) 
and Qo = 0.2, (—0.5, —0.5, 1), respectively. The dashed lines represent the approximate solution. The upper dotted line 
represents the value of the density contrast at turnaround for a spherical perturbation. 

Fig. 9, is the density contrast at virialization, while Fig. 10 represents the turnaround epoch for a prolate spheroid. 
The solid, short-dashed and long-dashed lines, represent respectively the time of turnaround for an isolated spheroid and 
for spheroids having Qo = 0.1, bi = (-0.5,-0.5,1) and Qo = 0.2, (-0.5,-0.5,1). The upper dotted line represents the 
value of the density contrast at turnaround for a spherical perturbation. 

2. 7. Collapse time and the MF 

In the previous subsections, we have studied the evolution of different mass elements and the calculation of collapse time. 

In fact, the problem of finding an expression for the MF is strictly connected to that of finding realistic estimates of 
collapse times of generic mass elements. Knowing this last, by means of statistical methods (as we shall see later) it is 
possible to obtain an expression for the MF. In the case of the spherical model (SPH), Zel'dovich approximation and 
ellipsoidal model, collapse times can be calculated analytically: in the case of spherical collapse (SPH) is simply 1.69/5;, in 
the case of Ze'ldovich approximation (ZEL) is I/A3, while the ellipsoidal collapse (ELL) can be calculated by finding the 
smallest positive root of Eq. (3.24) of M98 and by correcting for quasi spherical collapses as in Eq. (3.27) of M98 or using 
the analytical model of section 2.6. Second (2ND) and third (3ND) order ^ collapse times have been found by looking for 
the first instant at which J < 0, then using conventional root-finding algorithms to find the accurate value. 

Collapse times present the disadvantage that the passage from collapse to non-collapse takes place through infinity. 
It is more convenient to define the quantity F(q; A) as the inverse collapse time of the point q, when the initial field is 
smoothed at a resolution (mass variance) A: 

F is the basic dynamical quantity needed to construct the MF. In the SPH case F is simply proportional to the linear 
density contrast, F = Si/ 1.69; in the Zel'dovich (ZEL) case F is just equal to Ai. 

The quantity F is obviously a non-Gaussian process, and it is not a simple non-linear function of a Gaussian field 
(such as, for instance, a lognormal distribution): it is a complicated non- linear and non-local functional of the whole initial 
Gaussian perturbation field. Its 1-point PDF is the minimal amount of statistical information needed to construct the MF. 



^ The Zel'dovich approximation is the first term of a perturbative series; the perturbed quantity is not the density, as in Eulerian 

perturbation theory, but the displacement of the particles from the initial position. The problem of the evolution of a sclf-gravitating 
fluid can be reformulated in terms of equations for the displacement field. The equations for the displacement field have been found 
by several authors: Buchert (1989), Bouchet et al. (1995) and Catelan (1995). The Lagrangian system, (see Eqs. 1.42 and 1.43 in 
M98), can be perturbatively solved for small displacements. This has been done by Buchert (1989), Moutaxde et al. (1991), Bouchet 
et al. (1992), Buchert (1992), Buchert & Ehlcrs (1993), Lachiozo-Rcy (1993a,b), Buchert (1994), Bouchet et al. (1995) and Catelan 
(1995); see also Bouchet (1996) and Buchert (1996) for reviews. The first order solution, for suitable initial conditions (as given by 
the linear growing mode), is the well-known Zel'dovich (1970) approximation. The perturbative series has been calculated up to 
third order (from this comes the notation 3RD). 
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Fig. 11. Cumulative and differential PDFs of F. Figure taken from M98. 



The 1-point PDFs for the F processes, relative to different dynamical predictions, can be estimated by means of the 
Monte Carlo realizations (see M98). Fig. 11 shows such PDFs for two different spectral indexes, namely n = —2 and 1. 
Both cumulative and differential curves (the lattcrs in logarithm scale) arc shown. 

In the following, the ELL and 3RD predictions will be considered, and the PDFs will be mediated over four the spectral 
indexes, namely n = —2, —1, and 1. 

Summarizing, the MF problem can be formulated in terms of a non-Gaussian process F, representing the inverse 
collapse time of a generic mass element. If linear theory with a threshold (or spherical collapse) is considered, F is 
proportional to the linear density contrast, and is then Gaussian. More realistic estimates of collapse times lead to F 
processes which have a different statistical behavior, even in the rare event tail. 



3. Mass function. 

3.1. PS MF 

The mass function or multiplicity function can be described by the relation: 

dN = n{M)dM (92) 

that is the number of objects per unit volume, having a mass in the range M ed M+dM. The multiplicity function can also 
be used to define the luminosity function after having fixed the ratio ^ . Obtaining the mass function starting from that of 
luminosity is complicated since the ratio ^ is known with noteworthy uncertainty and it is different for different objects 
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Fig. 12. Comparison of PS multiplicity function with simulations. In the plot the solid line represents the multiplicity 
function obtained by PS, solid line, and the numerical multiplicity function obtained by YNY. Figure taken from Del 
Popolo (2005). 

and moreover the luminosity function for objects like galaxies depend on the morphological type (Binggeli & Tamman 
1985). Finally trying to determine the luminosity function observatively is problematic (see for example G.Efstathiou, 
R.S.Ellis, B.A.Peterson 1988). 

For the above reasons, the theoretical determination of the mass function is very important. 

The theoretical derivation of the mass function of gravitationally bound objects has been pioneered by Press & Schechter 
(1974). To incorporate the dynamics of the structure, PS formalism adopted top-hat spherical model, according to which 
the collapse condition for forming massive objects is determined purely by its local average density. To make statistical 
predictions, however, the PS formalism assumes that initial density field is Gaussian and that massive objects form in the 
peaks of the density field. 

This theory is based upon these hypotheses: 

— The linear density field is described by a stochastic Gaussian field. The statistics of the matter distribution is Gaussian. 

— The evolution of density perturbations is that described by the linear theory. Structures form in those regions where 
the overdensity linearly evolved and filtered with a top-hat filter exceeds a threshold 5c ((^c = 1-68, obtained from the 

spherical collapse model (Gunn & Gott 1972)). 

— for 6 > 6c regions collapse to points. The probability that an object forms at a certain point is proportional to the 
probability that the point is in a region with 5 > 5c given by: 




(93) 



The multiplicity function is given by: 
p{M,z) = -Po^dM = n{M)MdM 
If we add the conditions fl = 1, \5kf oc k'", the Press-Schechter solution is autosimilar and has the form: 



(94) 
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1 M , , :i+3l dM 



(95) 



where M*(2;) oc (1 + 2) "+3 . Several are the problems of the theory: 

— Statistical problems: in the limit of vanishing smoothing radii, or of infinite variance, the fraction of collapsed 
mass, asymptotes to 1/2. This is a signature of linear theory: only initially overdense regions, which constitute half 
of the mass, are able to collapse. Nonetheless, underdense regions can be included in larger overdense ones, or, more 
generally, non-collapsed regions have a finite probability of being included in larger collapsed ones; this is commonly 
called cloud-in- cloud problem,. PS argued that the missing mass would accrete on the formed structures, doubling their 
mass without changing the shape of the MF; however, they did not give a true demonstration of that. Then, they 
multiplied their MF by a "fudge factor" 2. Other authors (see Lucchin 1989) used to multiply the MF by a factor 
(1 + /), with / denoting the fraction of mass accreted by the already formed structures. 

— Dynamical problems: the heuristic derivation of the PS MF bypasses all the complications related to the highly non- 
linear dynamics of gravitational collapse. Spherical collapse helps in determining the Sc parameter and in identifying 
collapsed structures with virialized halos. However, the PS procedure completely ignores important dynamical elements, 
such as the role of tides and the transient filamentary geometry of collapsed structures. Moreover, supposing that 
every structure virializes just after collapse is a crude simplification: when a region collapses, all its substructure is 
supposed by PS to be erased at once, while in realistic cases the erasure of substructures is connected to the two-body 
interaction of already collapsed clumps, an important piece of gravitational dynamics which is completely missed by 
the PS procedure. 

— Geometrical problems: to estimate the mass function from the fraction of collapsed mass at a given scale it is 
necessary to relate the mass of the formed structure to the resolution 

In practice, the true geometry of the collapsed regions in the Lagrangian space (i.e. as mapped in the initial configura- 
tion) can be quite complex, especially at intermediate and small masses; in this case a different and more sophisticate 
mass assignment ought to be developed, so that geometry is taken into account. For instance, if structures are supposed 
to form in the peaks of the initial field, a different and more geometrical way to count collapsed structures could be 
based on peak abundances. 

Improvements of the model are due to several authors. Heavens and Peacock (1989), Peacock & Heavens (1990), BCEK, 
strictly connected to the lacking factor of 2 in the theory, which has been shown to be related to the so called " cloud in 
cloud" problem. Heavens and Peacock (1989) took into considerations the underdense regions {S < Sc)- Filtering a density 
field, S{x) on a given scale Rf, one obtains a set of points exceding the threshold S^- this set is named "excursion set" 
(see next subsection) of the filtered field (Adler 1981). In this field is possible to identify objects using special criteria (see 
Apple & Jones 1990). An object who had a given value S > Sc, at a fixed Rf, shall have a smaller value S = Sc, for a 
larger Rf. If Rf is furtherly increased it shall disappear under the threshold. As a consequence, objects belonging to the 
smaller scale of clustering shall be englobed in the larger one and as a consequence half of the mass is accounted. Heavens 
and Peacock (1989) solved the problem taking into consideration also the underdense regions with S < Sc- They used, 
differently from PS, the relation: 



where pc is the Gaussian distribution used by PS. This relation divides the points into two classes: those going from 
S > Sc to S = Sc, when Rf is increased, and they are associated to objects of radius > Rf, and points under the threshold 
having the probability Pup that in a following filtering, they have S > Sc- Using the relation for p^p obtained by Heavens 
and Peacock (1989) one gets: 



which solves the problem of the fudge factor of 2 in PS. 

In next subsections a wider description of the improvements to the PS theory are summarized. 

Despite all of its problems, the PS procedure for a long time proved successful, as compared to N-body simulations, 
and a good starting point for all the subsequent works on the subject (Efstathiou et al. 1988; Efstathiou & Rees 1988, 




(96) 



p{> Rf) = 2pG 



(97) 
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Fig. 13. The absorbing barrier problem. Figure taken by M98. 

BCEK), White, Efstathiou & Frank 1993, Jain & Bertschinger 1994, Lacey & Cole 1994, Efstathiou 1995, Bond & Myers 
1996b, Governato et al. 1999). 

Although the analytical framework of the PS model has been greatly refined and extended (see next subsection), 
more recently, it has been shown that the PS mass function, while qualitatively correct, disagrees with the results of high 
resolution N-body simulations. In particular, the PS formula overestimates the abundance of haloes near the characteristic 
mass and underestimates the abundance in the high mass tail (Efstathiou et al. 1988; White, Efstathiou & Frenk 1993; 
Lacey & Cole 1994; Tozzi & Governato 1998; Gross et al. 1998; Governato et al. 1999; YNY). The quoted discrepancy is 
not surprising since the PS model, as any other analytical model, should make several assumptions to get simple analytical 
predictions. In Fig. 12, it is plotted a comparison of PS multiplicity function with simulations. In the plot the solid line 
represents the multiplicity function obtained by PS, solid line, and the numerical multiplicity function obtained by YNY. 
The plot shows what is well known, namely: the PS mass function, while qualitatively correct, disagrees with the results of 
N-body simulations. In particular, the PS formula overestimates the abundance of haloes near the characteristic mass M* 
and underestimates the abundance in the high mass tail (Efstathiou et al. 1988; Lacey & Cole 1994; Tozzi & Governato 
1998; Gross et al. 1998; Governato et al. 1999). 

3.1.1. The cloud- in-cloud problem 

The cloud-in-cloud problem has origin in the following inconsistency of the original PS procedure. A collapse prediction is 
given to any point (of the Lagrangian space) for any resolution A; in other words, a whole trajectory in the Si A plane is 
given to any point, as in Fig. 13. Such trajectories start from at A=0 (vanishing resolution implies complete smoothing, 
and then vanishing density contrast), then wander around the plane, eventually upcrossing or even downcrossing the 
threshold line Si = Sc- When a trajectory lies above the threshold, the point is assumed to be part of a collapsed region of 
radius R' > R{A); it is clear that the size R'{A') of the formed structure is connected to the point A' of first upcrossing of 
the trajectory. On the other hand, when a trajectory downcrosses the barrier at a resolution A", the point is interpreted 
as not included in any region of size > R"{A") {R decreases with increasing A), which is clearly in contradiction with 
what stated above. 

To solve the cloud-in-cloud problem, a point whose trajectory has experienced its first upcrossing of the threshold line 
has to be considered as collapsed at that scale, regardless of its subsequent downcrossings. This can be done as follows: 
an absorbing barrier is put in correspondence with the threshold line, so as to eliminate any downcrossing event (BCEK). 
Alternatively, a non-collapse condition can be formulated as follows: a point is not collapsed at A if its density contrast at 
any A' < A is below the threshold (Peacock & Heavens 1990). 

The mathematical nature of the problem, and the resulting MF, strongly depend on the shape of the filter. For general 
filters, trajectories are strongly correlated in A, and then all the N-point correlations of the process at different resolutions 
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have to be known to solve the problem. However, if the smoothing filter sharply cuts the density field in the Fourier space, 
then independent modes are added as the resolution changes, and the resulting trajectories are Gaussian random walks. 
Such kind of filter is commonly called sharp k-space filter, it will be referred to as SKS filter throughout the text. 

In the SKS case, the problem is suitably solved within the diffusion framework proposed by BCEK (see next subsection). 



3.2. Excursion set approach 

In this section, we review the formulation of the excursion set approach, mostly following the treatment given by BCEK. 
The terminology "excursion set approach" was introduced by BCEK, to indicate that the MP determination is based on 
the statistics of those regions in which the linear density contrast 6i is larger than a threshold Sc (such regions are called 
excursion sets in the theory of stochastic processes; see, e.g., Adler 1981). The PS procedure is clearly included in this 
approach. This section presents those works which are based on the excursion set approach. 



3.2.1. Langevin equation 

As previously described, the statistical properties of the Gaussian field S{x) are completely specified by the two-point 
function in Pourier space, which is related to the power-spectrum P{k) by ((5(ki) 5(k2)) — {2n)^S^ (ki -I- k2)P(fci), where 
(5^ represents the Dirac delta function, and the brackets (•) denote ensemble averaging. Our Fourier transform convention 
is 5(k) = /dx(5(x)e*'^ ''. 

We want now to study the statistical properties of the density fluctuation field at some resolution scale Rf. This is 
introduced by convolving S{x) by some filter function T4^(|x' — x|, Rf), 

S{x, Rf) = J dx' W{\x' - x|, i?/) (5(x') = J dkW{kRf) S{k) e-'^'' , (98) 

where W is the Pourier transform of the filter. At each point x the smoothed field represents the weighted average 
of (5(x) over a spherical region of characteristic dimension Rf centred in x. The detailed properties of S{x,Rf) clearly 
depend upon the specific choice of window function. The most commonly used smoothing kernels are the top-hat filter 
Wth{\'^\, Rf) = 3 6(i?/ — |x|)/47ri?j, where 6(a;) is the Heaviside step function, and the Gaussian one WciXjRf) = 
(27ri?^)~^/^ exp(— a;^/2i?j). Recently, for convenience of analysis, top-hat filtering has been also applied in momentum 
space WsKs{k, Rf) — 6(fc/ — k), where kf = and kf — |k/|. This kernel is generally called sharp /c-space filter. 

While it is easy to associate a mass to real space top-hat filtering MxHiRf) — 47rpbi?j/3, there is always a bit of 
arbitrariness in assigning a mass to the other window functions. The most common procedure is to multiply the average 
density by the volume enclosed by the filter, obtaining Mc{Rf) = (27r)'^/^p{,i?^ and MsKs{Rf) = GTr'^PbkJ^ (Lacey & Cole 
1993 (LC93)). An alternative procedure, originally introduced by Bardeen et al. (1986) (BBKS), corresponds to the choice 
MsKs{Rf) — 47rphi?ip^/3, where Rth is chosen by requiring crsKsi^f) ~ ^th{Rth), and similarly for the Gaussian 
filter. In this way one obtains good agreement with numerical simulations of clustering growth (Lacey & Cole 1994). 

In order to mimic the accretion of matter one consideres a full hierarchy of decreasing resolution scales Rf (Peacock 
& Heavens 1990, Cole 1991, BCEK and LC93). The effect of varying Rf can be obtained by differentiating eq. 

This has the form of a Langevin equation, which shows how an infinitesimal change of the resolution scale Rf affects the 
value of the density fluctuation fleld (5(x, i?/) in the given position x through the action of the stochastic force ri{x.,Rf). 
In the limit Rf oo one has S{:x.;Rf) — s- 0, which can be adopted as initial condition for our flrst-order stochastic 
differential equation. Thus, by solving it, we can associate to each point x a trajectory S{x,Rf) obtained by varying 
the resolution scale Rf. Trajectories associated to different neighbouring points will be statistically influenced by the 
correlation properties of the force ry(x, Rf), i.e. of the underlying Gaussian field (5(x). On the other hand the coherence of 
each trajectory along the Rf direction depends exclusively on the analytic form of the filter function and vanishes for the 
sharp fc-space window (BCEK). With such a filter, by decreasing the smoothing length one adds up a new set of Pourier 
modes of the unsmoothcd distribution to (5(x, i?/). For a Gaussian field this is completely independent of the previous 
increments, and each trajectory S{x, Rf) becomes a Brownian random walk. 
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In the case of sharp fc-space fihering the notation greatly simphfies if we use as time variable the variance of the filtered 
density field, A = a^{kf) = ((5(fc/)2) = {2tt^)-'^ J^^' dkPP{k). In such a case the stochastic process reduces to a Wiener 
one, namely 



a(5(x. A) 



dA 

with (C(x,A)) 



C(x,A) , 
and 



(C(x,Ai)C(x,A2)) -<5,(Ai-A2) 



(100) 



(101) 



[see eq. (17) for the spatial correlation]. In the following we will adopt A as time variable, unless explicitly stated. The 
solution of the Langevin equation H100() in an arbitrary point of space (the position x is here understood) , with the initial 
condition 6{A = 0) = 0, is simply (5(A) — /J^ (iA'^(A'). By ensemble averaging this expression one obtains (S{A)) = and 
((5(Ai)(S(A2)) = min(Ai, A2), which uniquely determine the Gaussian distribution of S{A). 

As already told PS model is intrinsically flawed by the cloud-in-cloud problem, namely the fact that a fiuctuation on 
a given scale can contain substructures of smaller scales and the same fluid elements can be assigned, according to the PS 
collapse criterion, to haloes of different mass. Moreover, in a hierarchical scenario, one expects to find all the mass collapsed 
in objects of some scale, while the PS model can account only for half of it: this problem is intimately related to the fact 
that in a Gaussian field only half volume is overdense. Press and Schechter faced the problem by simply multiplying their 
result by a fudge factor of 2. In this section we review how the Langevin equations introduced above can be used to extend 
the PS theory in such a way to solve both problems. 

The solution of the cloud-in-cloud problem has been given by Peacock & Heavens (1990), Cole (1991) and BCEK. 
Their approach consists in considering at any given point the trajectory d{Rf) as a function of the filtering radius, and 
then determining the largest Rf at which S{Rf) upcrosses the threshold tf{zf) corresponding to the formation redshift Zf. 
This determines the largest mass collapsed at that point, all sub-structures having been erased. So, the computation of 
the mass function is equivalent to calculating the fraction of trajectories that first upcross the threshold t/ as the scale M 
decreases. The solution of the problem is enormously simplified for Brownian trajectories, that is for sharp /c-space filtered 
density fields. In such a case one only has to solve the Fokker-Planck equation for the probability density W{6, A) dS that 
the stochastic process at A assumes a value in the interval S,S + d6, 



dW{S,A) 



1 d^W{S, A) 
2 



dA 2 dS^ 
with the absorbing boundary condition yV{tf,A) 
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5d(5). The solution is well known 

(103) 



Defining S{A,tf) — J^^^d6W{S,A,tf) as survival probability of the trajectories, one obtains the density probability 
distribution of first-crossing variances by differentiation 



ldW{S,A,tf) 
2 dS 



t 



f 



V2nA^ 




(104) 



The function Vi (A) dA yields the probability that a realization of the random walk is absorbed by the barrier in the 
interval (A, A + dA) or, by the ergodic theorem, the probability that a fiuid element belongs to a structure with mass in 
the range [M(A -I- dA), M{A)]. Finally, the comoving number density of haloes with mass in the interval [M,M + dM] 
collapsed at redshift z/ is 



n{M,Zf)dM = ^Vi{A) 



dA 



dM 



dM 
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Inserting the expression of Vi (A) of eq. I|1U4|) in the latter equation one obtains the well-known PS expression for the mass 
function 

\2 • 
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Fig. 14. Examples of trajectories obtained by smoothing the li^p^^r density fluctuation field with a series of sharp fc-space 
filters (with decreasing resolution scales Rj = 1/fc/) at two points separated by the distance r. Here we consider r = 4Mpc 
and a CDM power spectrum (with density parameter = 1 and present Hubble constant Ho = 50kms~^ Mpc~^) linearly 
extrapolated to erg = 1. The heavy continuous line represents the trajectory associated to a point x in Lagrangian space. 
The trajectory associated to the point y such that |y — x| = r is plotted with a light continuous line. The dotted line is 
obtained by: i) considering the trajectory associated to y, ii) artificially removing any correlation with the trajectory at x 
and Hi) suitably rescaling the result. In practice, the dotted line represents a trajectory which is completely independent of 
the one associated to x, but which has the same statistical one-point properties. This is what is generally used in building 
up merger trees of dark matter haloes. 



The original fudge factor of 2 of the PS approach is now naturally justified. 

Previous investigations (e.g. Peacock & Heavens 1990; BCEK) have shown that only for sharp fc-space filtering it is 
possible to write an analytic formula for the mass function obtained from the excursion set approach. Numerical solutions 

of the cloud-in-cloud problem with physically more acceptable smoothing kernels like Gaussian and top-hat result in mass 
functions that are a factor of two lower in the high-mass tail and have different small-mass slopes compared with the PS 
formula. The standard interpretation of this result is that the excursion set method is not reliable for M <C M*, where 
M*, defined by A(M*) =t'j, is the typical mass collapsing at 2;/. 

3.2.2. Merging histories 

Merging histories are an important piece of information in the formation history of dark matter objects. They are a natural 
outcome of MF theories: a mass point, found in an object of mass M at a given time, will be found in another, more massive 
object at a subsequent time; from the conditional probabilities connected with such events it is possible to construct the 
statistics of accretion and merging histories of collapsed structures. This was attempted first by Carlberg (1990), whose 
results arc in contradiction with more recent works outlined in the following. Bower (1991) constructed merging histories 
by using the PS formalism, obtaining the same results as that obtained by means of the diffusion formalism, which are 
now reported. 

Suppose that we have smoothed the initial density distribution on a scale R using some spherically symmetric window 
function WM{r), where M{R) is the average mass contained within the window function. There are various possible 
choices for the form of the window function (c.f. LC93). We use a real-space top-hat window function, WM{r) = Q{R — 



28 



r)(47r_R'^/3) ^, where Q is the Heaviside step function. In this case M = inpoR^ /3, where po is the mean mass density of 
the universe. The mass variance S{M) = a'^{M) may be calculated from 

a^{M)^ ^ J P{k)W^{kR)k^dk, (107) 

where P{k) is the power spectrum of the matter density fluctuation, and W{kR) is the Fourier transform of the real space 
top- hat. 

The "excursion set" derivation due to BCEK leads naturally to the extended Press-Schechter formalism. The smoothed 
field S{M) is a Gaussian random variable with zero mean and variance S. The value of 5 executes a random walk as the 
smoothing scale is changed. Adopting an ansatz similar to that of the original Press-Schechter model, we associate the 
fraction of matter in collapsed objects in the mass interval M, M + dM at time t with the fraction of trajectories that make 
their first upcrossing through the threshold uj = dc{t) in the interval S, S + dS. This may be translated to a mass interval 
through equation H107() . The threshold 6c{t) corresponds to the critical density at which a pertubation will separate from 
the background expansion, turn around, and collapse. It is extrapolated using linear theory, dc{t) = Scfi/Diin(z), where 
-Diin(z) is the linear growth function. The halo mass function (here in the notation of LC93) is then: 
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The conditional mass function, the fraction of the trajectories in halos with mass Mi at zi that are in halos with mass 
Mo at zo (Ml < Mo, zo < zi) is 



/(S'i,cji I 5o,cjo)d5i = 
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The probability that a halo of mass AIq at redshift zq had a progenitor in the mass range (A/i , Af i 
(LC93): 
dP 
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dMi) is given by 
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where the factor M^/Mi converts the counting from mass weighting to number weighting. 

We can now derive two more useful quantities. Given the mass of a parent halo and the redshift step zq 
average number of progenitors with mass larger than M/ is: 



N = {Np{M \Mo)) = 
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We can also readily calculate the average fraction of Mq that dwelt in the form of progenitor halos of mass M > Mf. 
dP 



fp 



dM— (M, zi|Mo,zo). 
M, dM 
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and define the complimentary quantity for the average fraction of Mq that came from "accreted" mass, /acc — ^ — fp- 

In order to describe the accretion history of an object, it suffices to lower the barrier in a continuous way, and follow 
the position in A of the first upcrossing point: if this performs a discontinuous jump (which happens when the trajectory 
goes down and then up again), the object containing the mass point considered suffers a discontinuous merging with a 
structure of comparable size, while if the point moves continuously the object is just accreting material. It is clear that, if 
the trajectory is a random walk, the first upcrossing point will always perform discrete jumps; continuous accretion will 
be recognized only if an (arbitrary) minimum resolution step is fixed. 

LG93 also proposed a Monte Garlo approach to simulate ensembles of formation histories, based on realizing a large 
number of random walks. Their Monte Garlo method for simulating merging histories is commonly used to model the 
formation of virialized galactic halos, in which gas dynamical is inserted "by hand" . Such Monte Garlo models of galaxy 
formation will be discussed in §5.1.2. As a matter of fact, they found a weak inconsistency in their formalism (a probability 
density going slightly negative), probably caused by the simplicistic mass assignment. The same authors (Lacey & Gole 
1994) compared their results to N-body simulations, finding an overall satisfactory agreement. 
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Finally, Sheth (1995,1996) obtained a complete analytical description for the merging histories of objects formed from 
a Poisson distribution of seed masses, the problem analyzed by the original PS paper. The resulting MF has turned up to 
be the same as the distribution function proposed by Saslaw & Hamilton (1984). 

3.3. The statistics of the collapsed regions 

An idea which traces back to Doroshkevich (1970) is to suppose that structures form in the peaks of the initial density field. 
This idea became a standard paradigm in the framework of biased galaxy scenarios. Kaiser (1984) noted that high-level 
peaks show an enhanced correlation with respect to the underlying matter field, a fact which provided an explanation for 
the large correlation length of clusters with respect to galaxies, and gave freedom to tune the normalization of the CDM 
model to reproduce the large-scale distribution of galaxies. Peacock & Heavens (1985) and BBKS calculated a number of 
statistical expectation values for the peaks of a Gaussian random field, as the number density of peaks of given height. 
This quantity seemed suitable to determine a peak MF, but two important problems, recognized by BBKS hampered such 
a determination (see next subsection for details). 

In order to get the peak MF, it is necessary to study the stochastic properties of the collapsed regions, defined as 
regions where the linearly evolved density contrast 6{M, x) exceeds a threshold Sc- The correct choice of the threshold 6c, 
insofar as it is well defined, is a matter of debate. Many authors take the value 1.69 inspired by the spherical collapse 
model. On the other hand, comparison of the quasi-linear estimate of n(> M) (described below) with estimates from 
numerical simulations suggests a smaller value, Carlberg and Couchman (1989) advocating 6c = 1-4:4: and Efstathiou and 
Rees (1988) advocating 6c = 1.33 (but see also Brainerd and Villumsen (1992) and Katz, Quinn and Gelb (1992)). The 
collapsed regions occupy a volume fraction V given by the Gaussian distribution, 

f- = ^e-V2 (US) 
leading to 

V{iy)=eric{iy/V2)/2 (114) 
=(27r)-V2i.-ie--V2(i _ ^-2 ^ o(j,-4)) (115) 

To say more one needs to know the shape of the spectrum of perturbations. We shall list the relevant results given by 
BBKS. They involve only two moments of the spectrum, defined by 

{k''{M))=a-^{M) / k''exp{-k'^R})P{k)^ (116) 
Jo ^ 

{k\M))^a-\M) / k^eM-k^R})P{k)^ (117) 
Jo ^ 

The quantity (fc^) is the mean of the operator, ie., of the quantity S^^W^6. Similarly, (fc"') is the mean of V*. 

A relevant length scale is defined by i?J = 3{k^)/{k^). For any spectral index n > —1, it is easy to show that in the 
limit of small filtering scale Rf, 

Rf \l + nj ^ ^ 

For the case of CDM with .7 < n < 1, the ratio is in the range 1 to 3 for the entire range of cosmologically interesting 
masses. 

Another relevant length scale is (A:^)~^/^. On large filtering scales, such that P(fc)is increasing fairly strongly at 
k~^ ~ Rf, the ratio {k'^)~^^^ /Rf is close to 1. As the scale is reduced it increases, but is < 10 for M > IO^Mq. 
Finally, it is convenient to define the dimensionless parameter 

7(M) = {k')/{ky/' (119) 

It falls from about .7 to about .3 as M decreases from lO^^Mg to lO'^Mg, for .7 < n < 1. 

For sufficiently large u, each collapsed region is a sphere surrounding a single peak of 6. However, the departure from 
sphericity is considerable in the cosmologically interesting regime. BBKS show that a quantity x~^, which is roughly the 
fractional departure from sphericity, is well approximated by 

a; = 71/ + 6'(7, 7z/) (120) 
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where 

, ^ 3(l-7^) + (1.216-.97^)exp[-7/2(7./2)^] 
[3(1 - 72) + .45 + (7z//2)2] V2 + ^^/2 

We emphasize that this is the asphericity seen in the hnearly evolved, filtered density contrast. The asphericity in the 
true, unfiltered density contrast will be bigger, and will increase during collapse. 

Three useful number densities are given by BBKS. First, the density of the Euler number of the surfaces bounding 
the collapsed regions is 

Second, the number density of upcrossing points on these surfaces is 

mm"' 



(fc')) = ' i^' - 1)^-'^ (122) 



n^p{y,{k ),7) = — 

e-"'/2 (123) 



V — 1 



572(1-572/9)1/2 

An upcrossing point on a surface of constant 5 is defined as one where points along some arbitrarily chosen reference 
direction. 

The third number density is ripoak, the number density of peaks which are more than v standard deviations high. 
BBKS give expressions for ripeak, but they point out also that in the cosmologically interesting regime it is quite well 
approximated by riup. We shall use this approximation in what follows. It suggests that if a collapsed region contains 
several peaks, they are not buried deep inside it; rather, the boundary of the region is presumably corrugated, wrapping 
itself partially around each peak. 



3.3.1. The number density n{> M) 

The main application of these results is to estimate the number density n{> M) of gravitationally bound systems with 
mass bigger than M, at a given epoch before z^i{M) (where nl stands for "non linear"). The systems arc supposed to 
be identifiable by looking at the linearly evolved density contrast S{M,x). Each collapsed region, defined as one in which 
5(M, x) > 6c, is supposed to contain one or more systems with mass bigger than M. 

If each collapsed region is identified with a single system, then n(> M) = ricoii- In general this recipe is useless for 
lack of an expression for ricou- A different prescription, which does lead to a calculable expression, is to identify each peak 
within a collapsed region with a different collapsed object, 

n(> M) = npeak riup (124) 

This estimate (usually without the simplifying second equality) is widely used in the literature. It is certainly the same as 

the estimate n{> M) = ricoii for large v, where we know that there is just one peak per collapsed region. To what extent 
the prescriptions are the same for lower v is not known, because the number of peaks per collapsed region is not known. 

By means of the previous quoted theory it seemed suitable to determine a peak MF, but two important problems, 
recognized by BBKS, (who did not attempt to determine an MF from the peak number density) hampered such a 
determination: (i) the peak number density was based on the initial field smoothed at a single scale, so the peak MF 
suffered from the same cloud-in-cloud problem (the peak-in-peak problem) as the PS one; (ii) it was not clear which mass 
had to be assigned to a peak. 

The first (i) problem can be easily described as follows. If at some epoch the linearly evolved density contrast does 
have many peaks within a collapsed region, an interesting situation arises. At a somewhat earlier epoch, d{M,x) was 
smaller, and a separate contour d{M,x) = 5c was wrapped around each peak. In other words, each peak of the linearly 
evolved density contrast, filtered on scale M, was inside a single collapsed region, and presumably represented a separate 
gravitationally bound system. At the later epoch when the collapsed region encompasses many peaks of the linearly evolved 
density contrast, we have a bigger gravitationally bound system. If the original systems survive, the identification of each 
peak with a separate system is correct, but it misses the larger system which contains the original systems. Of course, 
missing this one system does not affect the total count much, so if this case is typical of collapsed regions containing 
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many peaks the estimate n{> M) — ripoak is better than the estimate n(> M) — ricou- If, on the other hand, the original 
systems have merged, that identification is wrong, and the whole of the collapsed region should be identified with just 
one gravitationally bound system. If this case is typical, the estimate n(> M) — ncoU would be better, if only we had a 
formula for it. Which case is the more likely? A clue is provided by the observation made earlier, that if there are several 
peaks in a collapsed region they typically seem to lie near the surface of the region, a part of the surface wrapping itself 
around each peak. This picture would suggest that the estimate n(> M) = ripeak is the more reasonable, the peaks of the 
linearly evolved density contrast in a typical collapsed region representing structures which have not existed long enough 
to merge. 

Even if the determination of a peak MF has some problems, several authors tried to formulate it in different ways, e.g. 
Colafrancesco et al. (1989). Let's Using the method of the last paper, let's denote by iVpk(j^, M)dv the number density of 
peaks of heightbetween v and v+dv and 

/"OO 

npk= / N^k{v,M)dv (125) 

In order to obtain the mass multiplicity for objects belonging to a catalog of contrast be-, one considers peaks on a 
hierarchy of resolution scales, M. Following the original PS choice, the mass of different objects is identified with the 
filtering mass M. The mass per unit volume n{M)MdAI , is obtained by differentiating with respect M the mass of the 
peaks of the filtered field with overdensity exceeding Sc, that is: 



n{M)MdM = (1 + /) 



d[np^{Sc;M)Mpi,{S„M)] 



dM 



dM. (126) 



which s a generalization of PS Ansatz, where 1 + / takes account of the secondary infall of matter initially in underdense 
regions, and A/pk is the average mass of peaks above 5c (see Eq. 3-7 of Colafrancesco et al. 1989) calculated modelling the 
peak as an ellipsoid, and estimated its mass by means of the volume inside the ellipsoid surface with density larger that 
a given threshold. 

As a matter of fact, there is not a general agreement on the actual validity of the peak paradigm. From the theoretical 
point of view, structures are not predicted to form in the peaks of the initial field; for instance, according to Zel'dovich 
approximation, structures form in the peaks of the largest eigenvalue of the deformation tensor. Then the peak paradigm 
can not be valid in general, except for the highest peaks. Some numerical simulations (Katz, Quinn & Gelb 1993; van de 
Weygaert & Babul 1994) have shown that galactic-size peaks often disrupt or merge with larger structures, as a result 
of tidal interactions with external structures. Manrique & Salvador-Sole have argued that such results are due to the 
lack of correction for the peak-in-peak problem. On the other hand. Bond & Myers (1996b) have found their peak-patch 
structures, which account for the peak-in-peak problem, to represent well N-body structures. 



The excursion set and peak approaches are somewhat complementary. In fact, excursion sets are effective in determining 
the total fraction of collapsed mass, and then the global normalization, but are not accurate in deciding how the collapsed 
mass fragments into clumps, i.e. to count the number of objects formed. On the contrary, the peak approach clearly 
determines the number of formed objects, but does not determine the mass to be associated with the structures, and 
hence the global normalization. 



3.3.2. Relations between the peak model and PS formalism 

The peak model can be connected to PS formalism. As known Press and Schechter (1974) worked with the differential 
number density, 

dn d , , , , 

dM - dM^O (12^) 

At a given epoch, if the filtering mass M is increased by an amount dM then v is increased by an amount dv, and the volume 
fraction occupied by the collapsed regions is reduced by an amount dV given by Eq. (|113|l . Press and Schechter suppose 
that the eliminated volume consists of objects with mass between M and M -fdAf, corresponding to the idealisation that 
filtering the density contrast on any mass scale M cuts out precisely those objects with mass less than M while leaving 
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unaffected objects with mass bigger than M. Ignoring the overdensity ; 
the number density dn of such objects is given by 



"^dM = 



M 



d{R}) 



dM 



du dV dn 



2Rj 



d{Rj) 



d((T2(M)) diy dV 
1 



(1 + 5c) of the collapsed regions this implies that 



(128) 



2(72 (M) 



2n 



1 

Vf 



V2 



(129) 



(130) 



Press and Schechter multiplied this formula by a factor 2, so that when integrated over all masses it would give the 
total mass density of the universe, rather than just the half corresponding to the regions of space where the linearly evolved 
density contrast is positive. They thus arrived at the estimate 



n(> M) ~ rip 



Jm 



2\t 



Gtt'^R 



-V e 
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dM' 



(131) 



In this equation, R'^ = Rf{M'), and similarly for (/c^)' and i^' . The factor 2 inserted by Press and Schechter is not justified 
(as previously described) by their argument, because the linearly evolved density contrast has nothing to do with reality in 
the non-linear regime a{M) > 1. On the other hand, the neglected overdensity gives a factor ~ (1 + ^c) = 2 to 3. Thus the 
factor 2 goes in the right direction, and the Press-Schechter formula is reasonably well founded theoretically. A somewhat 
different justification for the formula has been given by Bond et al (1991b). 

As told before, there are two alternative estimates n(> M) = Hps and n(> M) = npeak- A comparison of these two 
shows that the estimates agree to better than a factor 2 for < 2. Presumably, this indicates that in this regime the 
assumptions underlying the two estimates are compatible, in that increasing M by a small amount cuts out portions of 
the collapsed regions which have mass of order M and are centred on peaks with height of order i^(M). 

For large u, the Press-Schechter estimate falls below ripeak- This can be understood analytically, from the expression 



dn 



peak 



driup dn^p di/ dn^p d{k'^) duup dj 



dM dM dv dM a(fc2) dM d-i dM 
The first term dominates for large v, leading to the ratio 

/ \ 3/2 

dnps/dM _ / 3 \ _3 



dripeak/dM 



(132) 



(133) 



Apart from the factor 2, this is just the filter volume divided by the average volume of a collapsed region. 

As we saw earlier, the smallness of the size of the 'collapsed regions' is an artefact of the filtering. The conclusion is 
that for very rare fiuctuations, npeak is a better estimate that nps, the latter being considerably too small (Thomas & 
Couchman 1992). However, other sources of error are likely to be more important than the difference between Tipeak 

and 

nps- 



3.4. Dynamics and MF 

The models treated till now identifies collapsed structures as those regions whose linear density contrast exceeds some 
threshold. Many simplifications are used to obtain the mass function which at the end leads to neglects important elements 
such as the role of tides. Such simplifications could lead to oversimplified and misleading arguments. A number of authors 
have tried to insert elements of realistic dynamics in the theoretical MF. Such attempts are reviewed in the following. 

Some authors have inserted elements of realistic dynamics in the MF problem by extending the original PS approach 
or the difi^usion or the peak one (Lucchin & Matarrese 1988). This approach named PS-like method shall be analized in the 
next subsection. The peak-patch formalism by Bond & Myers (1996a), described above, is also characterized by a more 
realistic description of the dynamical evolution of peaks, even though structures are always identified through the peaks 
of the linear field. Other determinations of the MF were proposed by Henriksen & Lachieze-Rey (1990), where collapsed 
regions were identified by means of correlated velocity structures, and by Newman & Wasserman (1990) and Bernardeau 
& Schaeffer (1991), who related (in quite different ways) the MF to the correlation properties of the matter field. Monaco 
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(1995) constructed a MF in a PS-like approach, based on realistic collapse time estimates, found by means of extensions 
of the Zel'dovich approximation, and by the use of the homogeneous ellipsoid collapse model. 

One of the problems with the PS approach, reported by Cavaliere et al. (1991), is that it supposes matter clumps 

to instantaneously pass from non-collapsed to collapsed, and to be immediately incorporated in a larger clump. In other 
words, PS seems to imply vanishing time scales for the formation and destruction of clumps. As a matter of fact, PS 
simply does not contain any information on such timescales: the change of the MF with time is a combination of creation 
of new clumps, destruction of old clumps and accretion of mass onto existing clumps. Such terms cannot be disentangled 
by means of the PS approach alone, without further assumptions: for instance, the "static" procedure proposed by LC93, 
which is based only on statistics, cannot provide a precise definition of formation time (it is arbitrarily, though reasonably, 
defined as the time taken by a clump to double its mass). 

Cavaliere et al. (1991) proposed a dynamical procedure, based on creation and destruction time scales, to model the 
MF and an evolution equation for n{M,t). 

Blain & Longair (1993a,b) and Sasaki (1994) used a similar approach, obtaining identical results. 

A completely different approach to the MF was proposed by Silk (1978) and Silk & White (1978). Aggregation 
(and fragmentation) of collapsed clumps of similar size can be described by means of an aggregation kinetic equation 
(Smoluchowski 1916; Ernst 1986). The behavior of the MF given by Smoluchowsky equation has been reviewed elsewhere 
(see, e.g., Lucchin 1989; Cavaliere, Colafrancesco & Menci 1991b; Cavaliere, Menci & Tozzi 1994). 

More recently, Shaviv & Shaviv (1993, 1995) have analyzed the Smoluchowsky equation (always in the gravitational 
context) in a different way; at variance with Cavaliere and coworkers, they found their MF to depend on initial conditions. 
Another application of Smoluchowsky equation in a cosmological context is due to Edge et al. (1990), to explain the 
evolution of the X-ray luminosity function of galaxy clusters. 

As a general remark, such kinetic approaches can describe those aggregation events which take place between already 
collapsed clumps. The direct hierarchical (first) collapse of structures remains well described by a diffusion formalism, like 
the one proposed by BCEK. 

The first determination of a mass function, based on a self-consistent realistic dynamical approximation, is probably 
due to the works on the adhesion model. 

Some authors compared the predictions of clump formation, as given by the adhesion model, to N-body simulations, 
in ID (Doroshkevich & Kotok 1990; Williams et al. 1991) and 2D (Nusser & Dekel 1990; Kofman et al. 1992), finding 
satisfactory agreement. 

Vergassola et al. (1994), attempted an analytical estimate of the MF with adhesion, by using and extending a number 
of theorems demonstrated by Sinai (see references in the cited paper). They were able to find the asymptotic dependences 
of the MF: it behaves exactly like PS at large masses (the exact position of the typical mass is not determined), but has 
a different slope at small masses. However, they could not find the exact normalization of the MF. 

Lagrangian perturbations can be used to construct an MF fully based on realistic dynamics. At variance with adhesion 
theory, this dynamical approximation is of truncated type, i.e. small-scale power has to be filtered out to avoid small-scale 
multi-stream regions. Moreover, Lagrangian perturbations show interesting connections with the homogeneous ellipsoid 
collapse model, which can also be used to give reliable collapse time estimates, as in Monaco (1995). Such topics have 
been addressed in Monaco (1997a, b). 

3.4.1. PS-like approach 

In previous sections, we saw how finding realistic estimates of collapse times of generic mass elements. The problem of 
translating such information into an expression for the MF is of purely statistical nature. Such a statistical problem, in 
the simple case in which the inverse collapse time F is proportional to the initial density contrast, has received much 
attention in the scientific literature. Two main approaches were identified, namely the excursion set and the peak ones; 
the first approach was shown to be easier to manage than the second one, at the expense of a simplified treatment of the 
geometry of collapsed regions in Lagrangian space, while the peak approach, whose validity relies on the validity of the 
peak hypothesis, better takes into account geometry, at the expense of an increased complexity of the formalism, especially 
when trying to include a proper treatment of the peak-in-peak problem. 
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Fig. 15. PS-like n(A) curves for Ze'ldovich approximation and ansatz 1, compared to the PS one. Figure taken from M98. 
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Fig. 16. PS-like n(A) curves for ELL and ansatz 2, compared to the PS one and with a PS with (5c=1.5 and the fudge 
factor 2; see text. Figure taken from M98. 

In this section, the MF will bo determined by means of a simplified PS-likc, single-scale^ approach, which consists of 
estimating the probability of having initial conditions that make the mass element collapse. 

A first determination of the MF can be obtained by applying the same statistical approach as in the original PS paper: 
the fraction of collapsed mass is obtained by integrating, at a given fixed scale, over all initial conditions which make a 
mass element collapse before a given time. This determination obviously suffers from the same cloud-in-cloud problem 
as the PS one; however, most mass is predicted to finally collapse by realistic collapse time estimates, so a PS-like MF 
is nearly normalized, by more than 90%. It is then natural to suspect that an MF obtained by means of the absorbing 
barrier formalism can be not very different from the PS-like one, as only a minor part of the mass, lying in the strongest 
underdensities, has to be redistributed; this mass is not expected to influence the MF in any interesting mass range. As 
shown by M98, that under some conditions the PS-like and the absorbing barrier MFs are very similar; the simplified 
PS-like approach suffices in finding the main features of the dynamical MF (M98). 

Integration over initial conditions requires that initial conditions are specified, and that their joint PDF is known. In 
the PS case, where linear theory and spherical collapse were used, initial conditions were simply provided by the initial 
density contrast. In the most general case, initial conditions are given by the value of the density (or, equivalently, of the 
potential) at every point, so that a direct integration is hard to perform. An intermediate case is provided by the Zel'dovich 



35 



approximation, and by the other approximations which require the same initial conditions and eUipsoidal collapse. In this 
case, the joint PDF of initial conditions, the A eigenvalues, is known (Doroshkevich 1970): 

675V5 f 3 2 , 15 



Pa(Ai, A2, A3)dAidA2dA3 = ^^^3 exp (^-^Mi + 
X (Ai - A2)(Ai - A3)(A2 - \3)d\idX2dX3. (134) 

A is again the mass variance; /ii = Ai + A2 + A3 and fi2 — A1A2 + A1A3 + A2A3 are principal invariants of the Zel'dovich 
deformation tensor. It is convenient to express this PDF in terms of the linear density 61 and the x and y variables defined 
in Eq. H58|l: in this case the joint PDF is factorized into a Gaussian for Si and a joint PDF for x and y: 

P(Sux,y)dS,dxdy = ^=e^p^^—j dSi — ^ -j^e^p^~ — ix + xy + y ) 

X xy{x + y)dxdy = Ps,{Si; A)ddi x Px,y{x,y; A)dxdy. (135) 
The fraction of collapsed mass can then be obtained as follows: 



ni<A)^ dx dyP,Jx,y;A) dSi Ps.iSf, A), (136) 

Jo Jo Js^{x,y} 

where the function Sc{x,y), which substitutes the Sc parameter of PS, is defined as the solution of the equation: 

b,i6c,x,y)^bito), (137) 

where b{to) is the instant at which the MF is wanted (it will usually be set equal to one). By writing the function 6c as 
Sq — f{x,y), where Sq is the spherical value 1.69 and the positive function f{x,y), vanishing at the origin, gives the effect 
of the shear, it is possible to write the n{A) function as: 

n(A) = nps{A) x X(A), (138) 
where nps{A) is the PS curve, and 1(A) is a correction term: 

2^(A) = ^ / dx dyP^^y{x,y)eTipi-—f{x,y) + j-f{x,y) 



5q \ dx dy ^ 

n(A) curves have been calculated for Zel'dovich, ellipsoidal model, and for the two ansatze previously presented. Details 
of the calculations are reported by Monaco (1995). Fig. 15 presents the n(A) curve for Zel'dovich and for the first ansdtz, 
Eq. (|60|) . in which spherical collapse is recovered when Zel'dovich predicts a slower collapse. The canonical PS n(A)curve 
is shown for comparison; the PS curve has not been multiplied by the fudge factor of two at this stage, in order to compare 
the results of the PS-like procedures in the different cases, with no guarantee of normalization (in Monaco 1995 all the 
curves were multiplied by two). It can be seen that the Zel'dovich curve underestimates the number of large- mass'^ objects, 
and gives more intermediate- and small-mass objects, also thanks to its better normalization. The ansdtz curve reproduces 
the PS one at large masses, and reduces to the Zel'dovich one at small masses, as expected. 

Fig. 16 shows the ellipsoidal model prediction, in comparison with the second ansdtz, Eq. (|61|l (with e = 0.2), the 
canonical PS curve (without the fudge factor 2, as before) and a PS curve with 6c=l-5, representative of the typical 
outcome of N-body simulations (and then with the factor 2). Both the ansdtz curve and ellipsoidal model predict an 
overabundance of large-mass clumps with respect to the canonical PS curve: it is again demonstrated that a systematic 
displacement of collapse times from the spherical value influences the large-mass tail of the MF, even though spherical 
collapse is asymptotically recovered. In particular, the ellipsoidal model curve is quite similar to the PS 1.5 curve (PS 
curve with Sc = 1.5). It is however to be stressed that this similarity, while encouraging, has to be taken with care, as it 
is not clear how the collapsed objects predicted by this theory are related to the N-body clumps. As a technical remark, 
this ellipsoidal model curve has been computed by using the full be curves; in Monaco (1995) only the overdense curve 
was considered, and the Zel'dovich behavior was forced at large x and y values. 
^ I freely use the word mass in this context to indicate the large-mass (small A) or small-mass (large A) part of the MF 



36 



3.4.2. MF and n(A) 

The passage from the resolution to the mass variable requires knowledge of how collapsed matter gathers in clumps, an 
element which is missing in the excursion set approach; The volume of the excursion sets as a function of resolution, and 
then the mass of structures, is obtained by means of the rule: 



Mn{M)dM = Q 



dM 



dM = Qn{A) 



dA 



dM 



dM. (140) 
A more realistic resolution-mass relation would predict a whole distribution of masses to form at a given resolution: 
A^p(M;A). (141) 
The p{M; A) function gives the probability that a mass M is formed at a resolution A; its mean value will be of order: 

Mp{M; K)dM ~ QR{kf, (142) 



as the smoothing length -R(A) is the relevant scale in this case; the proportionality constant, of order one, will in general 
depend on the shape of the filter, on the power spectrum and on the resolution. The MF would then be given by: 

Mn{M)dM = p(yj^ n{A)p{M, A)dA^ dM , (143) 

i.e., by the n(A) curve convolved with the distribution p{M,A). 

Such a p distribution is expected, at small resolutions, to be peaked at its mean value, while its shape at large 
resolutions is expected be more complex, probably influenced by the details of the prescription chosen to fragment the 
collapsed medium (M98). Then, with respect to using the simple golden rule, the introduction of the p distribution is 
expected not to influence dramatically the large mass tail of the MF, while it is likely to influence the small-mass part, 
which is confirmed to be a not robust prediction of this kind of MF theories. Finally, differences in the p distributions for 
different filter shapes could in principle at least attenuate the differences between the MFs found with SKS and Gaussian 
filtering. 



3.4.3. Conclusions 

As shown in M98, a PS-like procedure is enough to obtain the main features of the MF. If SKS filtering is used, it has 
been demonstrated that the diffusion formalism can be extended to the non-Gaussian F process, by considering it as a 
diffusion process. The problem can then be transformed to the diffusion of a Wiener process with a moving absorbing 
barrier. For Gaussian smoothing, the simple approximation proposed by Peacock & Heavens (1991) has been shown very 
successful in finding the MF. 

The following conclusions can be drawn: 

1. A larger number of large-mass objects is expected to form with respect to the canonical PS prediction with (5c=1.69. 

2. The fact that spherical collapse is asymptotically recovered for the strongest overdensities does not guarantee that the 
"spherical" canonical PS MF, with (5c=1.69, is recovered at large masses. 

3. The large-mass part of the MF is considered robust with respect to the dynamical prediction. The ELL prediction 
tends to give more objects than the 3RD one in the large-mass tail; this is due to the fact that 3rd-order Lagrangian 
theory slightly underestimates spherical collapse; thus the ELL prediction is considered more believable in that range. 

4. PS-like gives fewer objects than the SKS one, by roughly a factor of 2, in the large- and intermediate-mass ranges. 

5. PS-like is very similar to the PS one with a, Sc — 1-5 parameter (indicated as PS 1.5). 

The small-mass part of the MF is not considered a robust prediction of the theory, for at least three reasons: 

1. The definition of collapse, which is based on the concept of orbit crossing, is not expected to reproduce common small- 
mass structures like virialized halos. OC regions rather represent those large-scale collapsed environments in which the 
virialized halos are embedded. 
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2. All the dynamical predictions used are considered good as long as the inverse collapse time is not small. Thus, the 
small-mass part of the MF is based on non-robust dynamical predictions. 

3. The p distribution of the forming masses, at a given resolution, is likely to significantly affect the small- mass part of 
the MF. 

The prediction of more large-mass objects, caused by the improved dynamical description of collapse, confirms some 
previous claims, for example by Lucchin & Matarrese (1988) and Porciani et al. (1996), who introduced non-Gaussianity 
of dynamical origin in the PS or diffusion approaches. This increase can be seen as the effect of tides on the dynamics of 
the mass element (Bertschinger & ,Jain 1994). Besides, even spherical collapse, when the global interpretation of collapse 
times is assumed (§2.3.4), leads to the prediction of more large-mass objects (Blanchard et al. 1992; Yano et al. 1996; 
Betancort-Rijo & Lopez-Corredoira 1996). 

The similarity of the (Gaussian-smoothed) dynamical MF with a PS one, with (5c=1.5, makes the dynamical MF be 
consistent with many numerical simulation (§2.2), even though this agreement is not a real proof of validity as (i) OC 
regions are not the halos extracted from simulations, and (ii) the resolution-mass relation is still treated in a simplified 
way. However, the dynamical MF theory does not predict the trend of lower 5c values at higher redshifts, observed in 
recent N-body simulations (see references in §2.2). An explanation of this behavior was proposed by Monaco (1995): the 
presence of small-scale power can effectively slow down gravitational collapse, through dynamical events of the kind of 
previrialization, proposed by Peebles (1990; see also Lokas et al. 1996; Bouchet 1996), or through dynamical friction with 
those possible particles which evaporate out of the collapsed structures (Antonuccio & Colafrancesco 1994). Such events 
could be effective for power spectra with n > — 1 (Lokas et al. 1996). In CDM-like spectra, the effective spectral index at 
M* is smaller than —1 at high redshifts, where (5c=1.5 (or even less) is found, but becomes larger at lower redshifts, where 
5c starts to increase to 1.7 or even larger values. However, such a behavior could also be due to a dependence of the p 
distribution on the power spectrum. 

3.5. More recent developments using the excursion set model: the moving barrier model and the 
multiplicity function 

As reported in the introduction, the PS model when compared to numerical simulations gives a smaller number of high- 
mass halos while giving a larger number of low-mass halos. The quoted discrepancy, that lead to research new expressions 
for the mass function, is not surprising since the PS model, as any other analytical model, should make several assumptions 
to get simple analytical predictions. The main assumptions that the PS model combines are the simple physics of the 
spherical collapse model with the assumption that the initial fluctuations were Gaussian and small. On average, initially 
denser regions collapse before less dense ones, which means that, at any given epoch, there is a critical density, 5c{z), 
which must be exceeded if collapse is to occur. In the spherical collapse model, this critical density does not depend on 
the mass of the collapsed object, while taking account of the effects of asphericity and tidal interaction with neighbors, 
it is possible to show that it is mass dependent (Del Popolo & Gambera 1998, SMT). In the second hand, the Gaussian 
nature of the fluctuation field means that a good approximation to the number density of bound objects that have mass 
m at time z is given by considering the barrier crossing statistics of many independent and uncorrelated random walks, 
where the barrier shape B{m,z), is connected to the collapse threshold. Simply changing the barrier shape, SMT showed 
that it is possible to incorporate the "quoted effects" ^ in the excursion set approach. Moreover, using the shape of the 
modified barrier in the excursion set approach, it is possible to obtain a good fit to the universal halo mass function 
^. The excursion set approach allows one to calculate good approximations to several important quantities, such as the 
"unconditional" and "conditional" mass functions. STl provided formulas to calculate these last quantities starting from 
the shape of the barrier. 

In the following, I'll use an improved version of the barrier obtained in Del Popolo & Gambera (1998) to get the 
multiplicity functions, which shall be compared with those obtained by PS, ST, JOl, YNY, and with numerical simulations 
of YNY. 

^ Namely that in the case of objects collapsing at the same time, the less massive regions must initially have been denser than 
the more massive ones. 

® Note that at present there is no good numerical test of analytic predictions for the low mass tail of the mass function. 



38 



In order to calculate the barrier shape, it is possible to follow Del Popolo & Gambera (1998) model. The equation 
governing the collapse of a density perturbation taking account angular momentum acquisition by protostructures can be 
obtained using a model due to Peebles (Peebles 1993) (see also Del Popolo & Gambera 1998, 1999). 

Let's consider an ensemble of gravitationally growing mass concentrations and suppose that the material in each system 
collects within the same potential well with inward pointing acceleration given by g{r) (see Del Popolo & Gambera 1998). 
We indicate with dP = f{L,rVr,t)dLdvrdr the probability that a particle can be found in the proper radius range r, 
r + dr, in the radial velocity range Vr — Vr + dvr and with angular momentum L = rvg in the range dL. The radial 
acceleration of the particle is: 

^.^-,(0 = ^-^ (144) 

Eq. (|144|l can be derived from a potential and then from Liouville's theorem it follows that the distribution function, /, 
satisfies the coUisionless Boltzmann equation: 



dt ^ dr dvr 



r3 



= (145) 



Assuming a non-zero cosmological constant Eq. (|144|) becomes: 
dvr GM L^ir) A ,^ 

(Peebles 1993; Bartlett & Silk 1993; Lahav 1991; Del Popolo & Gambera 1998, 1999). Integrating Eq. we have: 

1 fdrV GM f ^ A 2 ,,,,, 

- — h/ — 7r^dr+—r^ + e (147) 

2\dt J r J M^r^ 6 ^ ' 

where the value of the specific binding energy of the shell, e, can be obtained using the condition for turn-around, ^ = 0. 

In turn the binding energy of a growing mode solution is uniquely given by the linear overdensity, 5i, at time ti. From 
this overdensity, using the linear theory, we may obtain that of the turn-around epoch and then that of the collapse. We 
find the binding energy of the shell, C, using the relation between v and 5i for the growing mode (Peebles 1980) in Eq. 
(|147() and finally the linear overdensity at the time of collapse is given by: 



5c. = 5c, 



1+ / . +A 



&GM 



(148) 



where ai = 0.585, /3i = 0.46, a2 — 0.4 and /?2 — 0.02, where 5co = 1.68 is the critical threshold for a spherical model, 
Tj is the initial radius, rta is the turn-around radius, L the angular momentum. The angular momentum appearing in 
Eq. H148|) is the total angular momentum acquired by the proto-structure during evolution. In order to calculate L, it is 
possible to use the same model as described in Del Popolo & Gambera (1998, 1999) (more hints on the model and some 
of the model limits can be found in Del Popolo, Ercan & Gambera 2001). 

The CDM spectrum used to calculate the mass function plotted in Figs. 17-19 is that of BBKS (equation (G3)), with 
transfer function: 

= [11^(1 + 2-34g)] ^ 3 ^ 2 ^ (5 4g^)3 ^ (6.7i)4]-i/4 (149) 

2.34q 

(where q — fj^^r^j—^ ■ Here 9 = pa/ {l.68p.y) represents the ratio of the energy density in relativistic particles to that in 
photons {9 — 1 corresponds to photons and three flavors of relativistic neutrinos). The power spectrum was normalized to 
reproduce the observed abundance of rich cluster of galaxies (e.g., Bahcal & Fan 1998). 

In the excursion set approach, the average comoving number density of haloes of mass m the universal or "uncondi- 
tional" mass function, n(m,z), is given by: 

p dlogi^ 



dlogm 



(BCEK), where p is the background density, v = ^ t'tjri) ) ratio between the critical overdensity required for collapse 

in the spherical model, 6c{z), to the r.m.s. density fluctuation a{m), on the scale r of the initial size of the object m. The 
function vf{v) is obtained by computing the distribution of first crossings, f{v)di', of a barrier B{v), by independent, 
uncorrelated Brownian motion random walks. The mass function can be thus calculated once a shape for the barrier is 
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given and the power spectrum is known. In the case of spherical collapse, characterized by a constant barrier (for all v), 
PS and BCEK obtained: 



V \2 V 



(151) 



In the case of a nonspherical collapse, the shape of the barrier is no longer a constant and moreover it depends on mass 
(Del Popolo & Gambera 1998; SMT). As shown by STl, for a given barrier shape, B{S), where S = (^-^^ = ^ and 
(7, = \/S7 — 5 CO, the first crossing distribution is well approximated by: 



j(^s)ds = \ns)\cM-^)^ 

where T{S) is the sum of the first few terms in the Taylor expansion of B{S): 
(-S')"9"B(5) 



ns) = J2 



dS" 



(152) 



(153) 



The quantity Sf{S, t) is a function of the variable v alone, where v = {5c{t) / a(M))'^ . Since 5c and a evolve with time in the 
same way, the quantity Sf{S,t) is independent on time. Setting 2Sf{S,t) — vf{v), one obtains the so-called multiplicity 
function /(i^). That's why the shape of the barrier influences the form of the multiplicity function. 
In the case of the ellipsoidal barrier shape given in ST: 



B{a'^,z) = ^5c{z) 



{av) 



(154) 



the Eqs. (|152|l . (|153|l . give, after truncating the expansion at n = 5 (see ST): 

vf{v) = ^/av/2Ti[\ + (3{av'^)-°'g{a)] cxp {^Q.bav^l + I3{av^y°'f) (155) 
where 

-I 1 - « + - ... - I (156) 

Using the values for (3 and a of ST (a = 0.707, 5c{z) = 1.686(1 + z), /3 ~ 0.485 and a ~ 0.615) in Eq. one gets 

(STl): 

2 



v}{v) ~ Ai 1 + 



0.094 

{av) 



av 

' — exp{-aj/ 



1 



0.5 



{av)- 



/2} 



with j4i ~ 1. This last result is in good agreement with the fit of the simulated first crossing distribution (ST): 
vf{v)dv = A2 ( 1 + 

where p = 0.3, and a = 0.707. 

The normalizatition factor A2 has to satisfy the constraint: 



— exp{-av/2) 
Ztt 



f{v)dv = 1 

and as a consequence it is not an independent parameter, but is expressed in the form: 



(157) 



(158) 



A. = 



1 + 2-Ptt-^/^T{1/2^p) 



^ 0.3222^. 



(159) 



(160) 



If the barrier takes account of the cosmological constant, like in Eq. I|148() . using the same method that lead to Eq. 
H155|l . we have that: 

2 



vf{v) - ^1 1 + 



(3ig{ai) /?2.g(a2)\ [av 



{av)"- 



{av)' 



— exp {—av 
27r 



1 



/3i 



02 



{av) ^ {av) 



(161) 



In the case of the barrier given in Eq. (|148f) with A = 0, the "unconditional" multiplicity function can be approximated 



by: 

vf{v)^Ai (1 + 



{av) 



0.585 



av 

' — exp \—acv 

271 



{av) 



0.585 



(162) 
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where a = 0.707, b — 0.1218, c = 0.4019, d — 0.5526 and A4 ~ 1.75 is obtained from the normahzation condition. 

In the case of the barrier with non-zero cosmological constant, Eq. H148|l . a good approximation to the muhiphcity 
function is given by: 

2 



0.1218 0.0079 \ 



i^fiu) A 1 + — ™ + — ^ exp {-0.4019 



^ 0.5526 . 0.02 



[av) [av) 



} (163) 



where ^5 — 1.75. As previously reported, for matter of completeness, to the previous functions, namely PS, ST, Eq. (|163|l 
we have to add JOl, which satisfies the equation: 

vf{v) = 0.315ea;p(- | 0.61 + ln[a-^{M)] \^-^) (164) 

In order to express the above relation as a function of v, one substitutes a~^{M) = v/5c and I assume a constant value 
of 5c: that of the Einstein-de Sitter Universe namely 5c = 1.686. The above formula is valid for 0.5 < v < 4.8. 
YNY (Eq. 7, hereafter YNY7) proposed the following function to fit the numerical multiplicity function: 

vf{u) = Af,[l + {Bu/V2f]v^eM-{Byl^)\ (165) 
where, is a normalization factor to satisfy the unity constraint, f(v)dv = 1, therefore 

^6 - 2{B/^)^{T[D/2] + T[{C + D)/2]}-\ (166) 

The best-fit parameters are given as 5=0.893, C=1.39, and D=0.408, and from these parameters, Aq is constrained so 
that Ad = 0.298. 

This best-fit, function from Eq. 1)1651) . is shown in Figs. 17-19 and is only valid at 0.3 < < 3. 



3.5.1. Results 

The analytic multiplicity functions of PS, ST, JOl, YNY7, and Eq. (|163|l . are compared with the numerical simulations of 
YNY. Those simulations adopt the ACDM cosmological parameters of Vim = 0.3, Vt\ — 0.7, h — 0.7, and (Ts — 1-0, using 
512^ particles in common (see YNY for details). 

The numerical multiplicity functions are shown by crosses with errorbars in five panels of Fig. 17. All the data from 
the initial redshift to the present z = is compiled to draw the average curves (crosses) with error bars indicating the 
epoch to epoch variation. In the panel (f), all the numerical multiplicity functions are shown by thin lines. 

Four analytic multiplicity functions described in the previous section are also shown in this figure, that is PS (solid 
line), ST (dotted line), JOl (long-dashed line), YNY7 (dashed line), and the one of the Del Popolo (2005) described by 
Eq. (|163|) ( dot-dashed line). Since the data are available only in the region at < 3, these functions could be erroneous 
at > 3. 

Note that the comparison of the above curves, except for the PS model, with the results of N-body simulations show 
a very good agreement. However, there are some discrepancies between the YNY multiplicity function and other model 
functions (except the one in Del Popolo 2005). First, the multiplicity function of the Del Popolo (2005), similarly to that of 
YNY, in the low-i^ region oiv< 1, systematically falls below the ST and the JOl functions. In this region the multiplicity 
function of the Del Popolo (2005) is very close to that of YNY. 

As seen in Fig. 17, and in agreement with YNY, the numerical multiplicity functions reside between the ST and JOl 
multiplicity functions at 2 < < 3 (except for the run 35b). Additionally, the numerical multiplicity functions have an 
apparent peak a.t v ^ \ instead of the plateau that is seen in the JOl function. 

On the other hand, in the high-z/ region, where v is significantly larger than unity, the multiplicity function of the 
Del Popolo (2005) like YNY takes values between ST and JOl functions. These differences between numerical multiplicity 
functions and analytic ones, like ST, STl and JOl, are within 1 a error bars, and they are possibly due to the different box 
sizes adopted (see YNY for a discussion). To be more precise, throughout the peak range of 0.3 < < 3, the ST multiphcity 
function is in disagreement with the high mass resolution TV-body simulations of YNY and that of Del Popolo (2005). As 
shown by YNY the ST functional form provides a good fit to them only choosing parameter values of a = 0.664, p = 0.321, 
and A2 = 0.301. The multiplicity function obtained in the Del Popolo (2005) has a peak at '--^ 1 as in the ST function, 
and YNY numerical multiplicity function and YNY7, instead of a plateau as in the JOl function. 
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Ic 

Moreover, the functional form proposed in YNY, namely YNY7, provides a better fit when compared with the ST 
functional form but it is not based on theoretical background. The function obtained in this paper, similarly to YNY7 
provides a better fit to simulations than the ST functional form, and at the same time has been obtained from solid 
physical, theoretical, arguments. The better agreement observed between the multiplicity function of Del Popolo (2005) 
and YNY simulations, when compared with the ST, is connected to the shape of the barrier ((5c). As reported in Sec. 2, 
taking account of the effects of asphericity and tidal interaction with neighbors, Del Popolo & Gambera (1998), showed 
that the threshold is mass dependent, and in particular that of the set of objects that collapse at the same time, the less 
massive ones must initially have been denser than the more massive, since the less massive ones would have had to hold 
themselves together against stronger tidal forces. 
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Fig. 17. The multiplicity functions from five runs of YNY simulations (crosses with error bars) are shown in the panels 
(a-e), except for the panel (/) which shows the results from all the runs (thin lines). In each panel are the solid line 
represents PS multiplicity function, the dotted line the ST multiplicity function, the long-dashed line the JOl multiplicity 
function, the short-dashed line YNY7, and the dot-dashed line the one calculated by Del Popolo (2005). Figure taken from 
Del Popolo (2005). 

The shape of the barrier given in Eq. H148|l is a direct consequence of the angular momentum acquired by the proto- 
structure during evolution and the effects of the cosmological constant. 

Similarly to ST, the barrier increases with S (decrease with mass, M) differently from other models (see Monaco 1997a, 
b) . It is interesting to note that the increase of the barrier with S has several important consequences and these models 
have a richer structure than the constant barrier model. 

The decrease of the barrier with mass means that, in order to form structure, more massive peaks must cross a lower 
threshold, Sc{i^,z), with respect to under-dense ones. At the same time, since the probability to find high peaks is larger 
in more dense regions, this means that, statistically, in order to form structure, peaks in more dense regions may have a 
lower value of the threshold, Scii', z), with respect to those of under-dense regions. This is due to the fact that less massive 
objects are more influenced by external tides, and consequently they must be more overdense to collapse by a given time. 
In fact, the angular momentum acquired by a shell centred on a peak in the CDM density distribution is anti-correlated 
with density: high-density peaks acquire less angular momentum than low-density peaks (Hoffman 1986; Ryden 1988). A 
larger amount of angular momentum acquired by low-density peaks (with respect to the high-density ones) implies that 
these peaks can more easily resist gravitational collapse and consequently it is more difficult for them to form structure. 
Therefore, on small scales, where the shear is statistically greater, structures need, on average, a higher density contrast 
to collapse. 

It is evident that the effect of a non-zero cosmological constant adds to that of L. The effect of a non-zero cosmological 
constant is that of slightly changing the evolution of the multiplicity function with respect to open models with the same 
value of Qa. This is caused by the fact that in a fiat universe with f2A > 0, the density of the universe remains close to the 
critical value later in time, promoting perturbation growth at lower redshift. The evolution is more rapid for larger values 
(in absolute value) of the spectral index, n. 

As previously reported, the ST model gives a better fit to simulations than PS model, but it has some discrepancies 
with simulations. ST model was introduced at the beginning (Shcth & Tormen 1999) as a fit to the GIF simulations and 
in a subsequent paper (SMT) was recognized the importance of asphcrical collapse in the functional form of the mass 
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Fig. 18. The best-fit multiplicity function. In the plot the solid line represents the multiplicity function obtained in Del 

Popolo (2005), the short-dashed line YNY7, the dotted line the ST multiplicity function, the long-dashed line the JOl 
multiplicity function. The errorbars with open circles represents the run 140 of YNY, those with filled squares the case 
70b, those with open squares the case 70a, those with filled circles the case 35b, those with crosses the case 35a. Figure 
taken by Del Popolo (2005). 

function. The effects of asphericity were taken into account by changing the functional form of the critical overdensity 
(barrier) by means of a simple intuitive parameterization of elliptical collapse of isolated spheroids. The model proposed in 
Del Popolo (2005) has several similitudes with ST and STl models, namely it uses the excursion set approach as extended 
by STl to calculate the multiplicity function, but at the same time it differs from ST and STl for the way the barrier was 
calculated and for the fact that takes account of angular momentum acquisition, and a non-zero cosmological constant, 
things that are not taken into account into ST and STl. These differences gives rise to a multiplicity function in better 
agreement with simulations. This shows the importance of the form of the barrier. The improvement of the multiplicity 
function of Del Popolo (2005) and ST with respect to PS is probably connected also to the fact that incorporating the 
non-spherical collapse with increasing barrier in the excursion set approach results in a model in which fragmentation and 
mergers may occur, effects important in structure formation. 

In the case of non-spherical collapse with increasing barrier, a small fraction of the mass in the Universe remains 
unbound, while for the spherical dynamics, at the given time, all the mass is bound up in collapsed objects. Moreover, 
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Fig. 19. Time dependence of the multiplicity function from the 35a run, for four redshift ranges of < z < 1 (open 
circles), 1 < z < 3 (open squares), 3 < 2 < 6 (open triangles), and z > 6, (crosses). Also shown are YNY7 (solid line) and 
the model of Del Popolo (2005) (dot-dashed line). Figure taken by Del Popolo (2005). 



incorporating the non-spherical collapse with increasing barrier in the excursion set approach results in a model in which 
fragmentation and mergers may occur (ST). If the barrier decreases with S (Monaco 1997a, b), this implies that all walks 
are guaranteed to cross it and so there is no fragmentation associated with this barrier shape. 

In other words, the excursion set approach with a barrier taking account effects of physics of structure formation gives 
rise to good approximations to the numerical multiplicity function: the approximation goodness increases with a more 
improved form of the barrier (taking account more and more physical effects: angular momentum acquisition, non zero 
cosmological constant, etc). Another important aspect of the quoted method is its noteworthy versatility: for example it 
is very easy to take account of the presence of a non zero cosmological constant englobing it in the barrier. I recall that 
the YNY numerical multiplicity function assumes a non zero cosmological constant while the theoretical models (ST,ST1, 
JOl) does not take this into account. 

Fig. 19 shows the multiplicity function from the 35a run, for four redshift ranges of < 2: < 1 (open circles), 1 < z < 3 
(open squares), 3 < z < 6 (open triangles), and z > 6, (crosses). At high redshifts, high-i^ halos in the exponential part of 
the YNY7 (solid line) function and Eq. (|163|l (dot-dashed line of Del Popolo 2005) are probed. As redshift decreases, the 
probe window moves to the lower-;/ region. Fig. 19 shows that the multiplicity function of Del Popolo (2005), Eq. H163(l . 
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and YNY7 both gives a good fit to the numerical simulations. For small values of ly, Eq. (|163fl is a slightly better fit to 
data, and at large values of v the two functions decays in the same way. 

3.5.2 Multiplicity function evolution 

In this section, I compare the analytic mass function of the present paper with that of ST, and with ? (R03) simulation 
results at several redshifts. In Fig. 7, I compare the mass function of the present paper with ST, and with R03 simulation 
results at several redshifts. In the figure, the solid Hne represents the ST mass function at z = 0,5,8,15, going from 
right to left, respectively. The dashed line the mass function of the present paper for the same values of the redshift, the 
errorbars with open squares, crosses, open triangles and solid triangles represents R03 at z = 0, 5, 8, 15. Fig. 7 shows that 
the ST function provides a good fit to R03 data, except at very high redshifts, where it significantly overpredicts the halo 
abundance. At all redshifts up to z—10, the difference is ~ 10% for each of our well sampled mass bins. However, the 
ST function begins to overpredict the number of haloes increasingly with redshift for z>10, up to ~ 50% by z=15. The 
simulation mass functions appear to be generally steeper than the ST function, especially at high redshifts. This is in 
agreement with the theoretical mass function calculated in the present paper which gives a better description of the R03 
mass function for higher values of z for which the ST mass function overpredicts the simulation results. 

In Fig. 8, I plot the mass function for all of our outputs in the f{a) — ln(fT~-'^) plane. Large values of Incr^^ correspond 
to rare haloes of high redshift and/or high mass, while small values of Incr"^ describe haloes of low mass and redshift 
combinations. The solid line is the ST mass function while the dashed line the one obtained in the present paper and the 
dotted line represents a crude multiplicative factor to the ST function as follows, with Sco = 1.686: 



/(a) = /(a;ST) 



exp[~0.7/{a[cosh{2a)f)] 



(167) 



valid over the range of -1.7 < Incr^^ < 0.9. The ST and the mass function of the present paper differs more in the high 
mass region, where the mass function of the present paper is steeper than ST and in better agreement with numerical 
simulations data than ST mass function. The ST function fits the simulated mass function to better than 10% over the 
range of -1.7 < In cr~^ < 0.5 while it appears to significantly overpredict haloes for Incr"^ > 0.5. The magnitude of the ST 
overprediction at high values of Incr"^ is consistent with being a function purely of Incr^^ rather than redshift, a natural 
consequence of the fact that the mass function is self similar in time (e.g. ?; ?; ?). The empirical adjustment to the ST 
mass function (Eq. (|167|l ). dotted line, describes much better numerical simulations data: for -1.7 < In cr^^ < 0.5, Eq. 
H167|l matches R03 data to better than 10% for well-sampled bins, while for 0.5 < lncr~^ < 0.9, where Poisson errors are 
larger, data is matched to roughly 20%. 

YNY and ? made a similar choice, namely they introduced an empirical mass function obtained from a fit to their 
simulations that gives a better fit to simulations than ST model. It is important to stress that even if the functional 
forms proposed in R03, YNY and ? provide a better fit to simulations when compared with the ST functional form, 
they are not based on theoretical background. The function obtained in the present paper, similarly, for example, to R03 
provides a better fit to simulations than the ST functional form, and at the same time has been obtained from solid 
physical, theoretical, arguments. The better agreement observed between the mass function of the present paper and R03 
simulations, when compared with the ST, is connected to the shape of the barrier ((5c). 

In other words, the excursion set approach with a barrier taking account effects of physics of structure formation gives 
rise to good approximations to the numerical multiplicity function: the approximation goodness increases with a more 
improved form of the barrier (taking account more and more physical effects: angular momentum acquisition, non zero 
cosmological constant, etc). Another important aspect of the quoted method is its noteworthy versatility: for example it 
is very easy to take account of the presence of a non zero cosmological constant englobing it in the barrier. I recall that 
the YNY numerical multiplicity function assumes a non zero cosmological constant while the theoretical models (ST, ?, 
?) does not take this into account. 

4. PROSPECTS AND CONCLUSIONS 



In the previous sections, we have seen that from the seminal paper of PS to last results, large progress has been made 
in the field of the cosmological mass function. If the PS model, explained quite well the N-body simulations of 80' and 
first part of 90', its limits have been shown by more recent simulations (e.g. ?; YNY). Its theoretical limits, for example 
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Fig. 20. Comparison of the mass function evolution calculated in the present paper with ST mass function and R03 
simulations. Solid curves are the Sheth & Tormen function at z=0, 5, 8, & 15 (from right to left). The dashed line the 
mass function of the present paper for the same values of the redshift, the errorbars with open squares, crosses, open 
triangles and solid triangles represents R03 result at the same redshift. 




Fig. 21. Mass function plotted in redshift independent form for all of R03 outputs: redshifts used are 0, 1., 2., 3., 4., 5., 
6.2, 7.8, 10., 12.1, 14.5. The solid line is ST prediction while the dashed and dotted line represent the result of the present 
paper and Eq. H167(l . respectively. 

the lack of the "fudge factor of 2" , has been solved in some papers (e.g. ?). We have reviewed the excursion set model, the 
relation between peak model and the MF, one approach (PS-like approach) that takes into account the dynamics in the 
MF theory. All the previous results, give a MF in better agreement with simulations but only in more recent years several 
works have shown how it is possible to obtain MF functions in very good agreement with simulations (?, ?, ST, ?, ?, ?, 
YNY). 

Even if the last models for the mass function gives much better results than the previous ones, it is necessary to recall 
which are the limits of the theory and the problems that must be attacked. 
The main problems are the following: 

— It is necessary to add more physics to the dynamical MF theory, to describe events such as the aggregation and 
fragmentation of already collapsed structures. In this way, the dynamics inside structures could be resolved; this is nec- 
essary to describe objects such as galaxies inside larger structures, as groups or clusters. The theories presented in the 
last subsection are somehow able to take account of some of the quoted effects, like fragmentation, angular momentum 
acquisition, presence of a non zero cosmological constant. Changing the shape of the barrier it is possible to take into 
account more physical effects to the dynamical MF theory. 

— It is necessary to take in full account the geometry of collapsed regions in Lagrangian space, going beyond the usual 
golden rule described in previous sections. In the realistic case in which the stochastic process on which the MF is based 
is non Gaussian, this investigation can be performed through Monte Carlo simulations of initial fields. 
— It is necessary to understand the role of filtering in the MF problem, as different kinds of filters (sharp fc-space, Gaussian) 
lead to different and not equivalent formulations of the problem. 
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— It is necessary to understand in detail what is the total mass of a structure, and to find a rigorous, well-posed and easily 
interpretable definition for it, suitable to help in the interpretation of observations. 

Another important point to remark is connected to N-body simulations. For example, as shown by YNY there is a 
discrepancy in the numerical multiplicity functions from various simulation runs. There arc three strategies to resolve this 
discrepancy. The first is to run simulations having still higher mass dynamic range free from the box size effect. The second 
is to increase the number of realizations, because there is a scatter from the runs using the same box size. The third is 
to run simulations whose box size is smaller than that of the present work, although it might sound contradictorily. Prom 
simulations with smaller box size, one obtain the information on the conditional nniltiplicity function which coincides with 
the unconditional multiplicity function at v « 1. Comparing the unconditional multiplicity function from simulations 
with a large box size and the conditional multiplicity function from those of a small box size will offer not only the clues 
to resolve the above mentioned discrepancies, but also insights into the mechanism how the PS ansatz works to reproduce 
the numerical multiplicity function. 
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